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HIGH-SPEED FLOW PAST WINGS 


By Helge Nprstrud 
Lockheed- Georgia Company 


SUMMARY 


The analytical solution to the transonic small-perturbation equation which describes 
steady compressible flow past finite wings at subsonic speeds can be expressed as a non- . 
linear integral equation with the perturbation velocity potential as the unknown function. 

This known formulation is substituted by a system of nonlinear algebraic equations to which 
various methods are applicable for its solution. Due to the presence of mathematical dis- 
continuities in the flow solutions, however, a main computational difficulty in the present 
study was to ensure uniqueness of the solutions when local velocities on the wing exceeded 
the speed of sound. For continuous solutions this was achieved by embedding the algebraic 
system in an one-parameter operator homotopy in order to apply the method of parametric 
differentiation. The solution to the initial system of equations appears then as a solution 
to a Cauchy problem where the initial condition is related to the accompanying incom- 
pressible flow solution. 

In using this technique, however, a continuous dependence of the solution development 
on the initial data is lost when the solution reaches the minimum bifurcation point. A steepest 
descent iteration technique has, therefore, been added to the computational scheme for the 
calculation of discontinuous flow solutions. Results for purely subsonic flows and supersonic 
flows with and without compression shocks are given and compared with other available 
theoretical solutions. 


INTRODUCTION 


The importance of having aircraft flying in the transonic speed range either in the 
cruise mode or a maneuvering mode has been expressed by both civilian and military 
operators. Although the financial benefits from a commercial transonic airplane seems at 
present to be somewhat questionable, no definite conclusion can be drawn before a 
thorough understanding has been reached of the various operational problems involved. 

The technical difficulties, however, associated with the development of such an aircraft 
are well known to both the experimentalists and the theoretician in the area of predicting 
the aerodynamic loads and the operating boundaries. 

Early attempts to theoretically analyze compressible flow past finite wings at subsonic 
speeds beyond an application of simple sweep or strip theory were based on the equivalence 
theorem which concerns itself with geometrically equivalent bodies (ref. 1). A solution was 
thus obtainable for small aspect ratio wings, but required the solution to the equivalent 
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nonlifting axisymmetric case. A more general three-dimensional treatment of the transonic 
flow problem, however, seems to have been first given by Alksne and Spreiter (ref. 2) who 
extended the concept of local linearization to nonlifting wings, but were limited to oncoming 
flows which deviated little from the sonic speed (see reference 3 for a more recent and 
detailed exposition). Similar restriction on the freestream condition also applies to the work 
of Burg (ref. 4) which is based on still another linearized form of the governing transonic 
small-perturbation equation. Nprstrud (ref. 5) adopted the integral equation approach to 
the solution of the same governing nonlinear equation, but results were confined to subcritical 
flows. The analytical difficulties associated with the three-dimensional transonic flow problem 
were clearly a matter which was inherited from the more familiar planar case and a solution 
to the spatial problem dependent to a large extent on the progress made in solving the problem 
in two dimensions, see e.g. reference 6 for a comprehensive literature review of the subject. 

The first successful evaluation of mixed subsonic/supersonic three-dimensional flow 
with embedded shocks were recently given by Bailey and Steger (ref. 7) who numerically 
solved the steady state small-perturbation equation using relaxation techniques. Their 
approach includes a hybrid combination of the perturbation velocity potential and the per- 
turbation velocities as the dependent variables and calculations for lifting, supercritical 
flows are presented. Similar relaxation schemes have also been applied by Caradonna and 
Isom (ref. 8) to the transonic flow problem of hovering helicopter blades at zero lift. Still 
further numerical investigations of the transonic wing alone or wing-body problem are found 
in references 9, 10 and 11 and this has brought the subject- of transonic flow theory to an 
initial stage of overall completeness. 

The present study refines and extends the work of reference 5 to lifting flows with the 
inclusion of shock discontinuities in three-dimensions. The approach taken follows some 
fundamental steps proposed for two-dimensional flows by Oswatitsch in 1950 (ref. 12) and it 
can be described as semi-analytical . The formal integral solution to the transonic small- 
perturbation equation with associated boundary conditions is regarded as more amendable to 
numerical analysis than the differential formulation. This is especially true for subcritical 
flows. For supercritical flows, however, some additional considerations must be given 
to guarantee the uniqueness of the flow solution. In either case, the governing integral 
equation is first replaced by a system of nonlinear algebraic equations and then alternative 
methods of solution to the algebraic system are given. 


SYMBOLS 


A. ,B.,C. 

i i i 

AR 

Ci 

c 

E(ix-?|,|Z-C|;r) 

F. 


velocity parameters 

aspect ratio 

cosine integral 

speed of sound 

influence function 

functional defined in equation (23) 
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f. 

I 


G. 

i 

H 

v 


J 

M 


m (x , Z) 
N v 


q(5,C) 

R 

r(x,Z) 

S 

Si 


s 

U,V,W 


u,v,w 

Y,Z 


x,y,z 


Greek Letters: 

a 

y(?,C) 

2 

v 

A 

6 .. 

1 1 

€ . . 
m 

H 


U 


P / a / X 


functional defined in equation (37) 

operator defined in equation (30) 

Struve function of integer order v 

Jacobian matrix 
local Mach number 
scaling parameter, see equation (22) 
Neumann function of integer order v 

source distribution, see equation (8) 

distance function 

velocity decay parameter 

wing planform area 

sine integral 

specific entropy 

transformed perturbation velocity components 
velocity components 

2 1/2 

transformed coord inates= (1 - M ) y,z 
Cartesian coordinates 

function describing wing surface distribution 


angle of attack 

vorticity distribution, see equation (9) 

Laplace operator 
incremental value 
Kronecker delta 

influence coefficients 

ratio of specific heats 
differentiation parameter 
integration variables 
arguments 
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ANALYTICAL ANALYSIS 


Governing Differential Equations 

Three-dimensional irrotational compressible flow of an ideal gas around a finite 
wing generates a velocity field (u,v,w) = grad I which is defined in a body-fixed Cartesian 
coordinate system x,y,z. With the further assumption of small disturbances in the plane 
normal to the direction of the oncoming steady flow (i .e . , in the yz-plane) the governing 
partial differential equation for the velocity potential 1= §(x,y,z) will be written as 


? .. 

(1 - M )$ + § + § = ‘ (h + 1 ) [ I - u ]$• 

» xx yy zz u x » xx 


0 ) 


where M (= u /c ) is the freestream Mach number and.x designates the ratio of specific 

00 CO GO 

heats. Equation (1) is a particular form of the set of equations which can be deduced from 
the transonic small-perturbation theory and it can be transformed to 


<P. 


xx 


Cpyy + ®P. 


ZZ 


CD CD 

Y X Y XX 


( 2 ) 


after a normalized velocity perturbation potential cp = cp (x,Y,Z) is defined as 

(x + 1) 

cp = = ($-u x-v y-w z) 

' / 00 00 ’ 00 

(1 -M> 


(3) 
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The employed coordinate transformation Y,Z = (1 - Mj ' y,z corresponds to a form of the 
Prandtl-Glauert transformation and we are thus limiting the present analysis to far-fields 
with subsonic flow (M < 1). This restriction is arbitrary since a different coordinate trans- 
formation would make a similar perturbation analysis valid for freestream Mach numbers 
M 2 1 . 

CD 

The introduction of a velocity potential $ in equation (1) is equivalent to imposing 
the condition of irrotationality 

rot(grad f) = 0 (4) 


on the flow field and, hence, equations (1) and (4) must be regarded as the system of equa- 
tions to be solved under given boundary conditions. Since equation (4) implies a charac- 
terization of the subject flow as global isentropic (s = const.) and assumes no vorticity in 
the upstream far field region, the entropy increase across any embedded shock surface is 
regarded negligible. Thus, the entropy condition which forbids the occurrence of expansion 
shocks must be incorporated in the analysis by means other than the use of the specific 
entropy measure . This missing formulation will be dealt with in the numerical analysis section 


Boundary conditions . - Let the planform projection of the wing be situated in the plane 
y = 0 and let the distribution y* = h ± (x,z) define the wing upper^ ^ and lower ^ ^ surface, 


w 


respectively. The simplified boundary condition expressing tangential flow at the wing 
surface is then 


y - ±0: $ “ u h (x,z) 

y =° X 


(5a) 


or in fransformed variables 

Y = ±0 : 

-1 


(k + 1)M' 


<Py 


(1 - M ) 


W 2 


[h (x, z) - tana] 


(5b) 


where a(= tan (v ^/u )) designates the angle of attack which will be limited to small 

values. Furthermore, equation (5b) assumes w^ = 0 for simplicity of the analysis . The 

condition at infinity shall be that of vanishing flow disturbances and will be stated as 


<x 2 + Y 2 + Z 2 ) ,/2 


cp = 0 


(6) 


In the case of circulatory (or lifting) flows, the Kutta condition at the trailing edge 
of the wing must be taken into consideration in order to ensure finite velocities of the 
flow as it leaves the wing surface. The associated potential jump Acpatthe trailing edge 
varies in a fashion similar to the spanwise load distribution, but remains constant in the 
streamwise direction in the so-called vortex wake behind the wing. 


The harmonic solution . - For the purpose of constructing a solution to the nonlinear 
partial differential equation (2) we will first consider the problem of obtaining a harmonic 
solution to the Laplace equation 

' p xx +T VY + 'f’zZ = 0 (7) 

This equation is well known from linearized subsonic theory and we can identify a solution 
cp = cp(x,Y,Z) which satisfies equation (7) and the boundary conditions (5b) and (6) as a 

2 2-1 

transformed Prandtl-Glauert potential with the conversion factor (h + 1)M (1 -M ) , 

CO CO ' 

see equation (3) . 

Analytical solutions to equation (7) for given wing configuration (or boundary 
conditions) are rare and, in fact, it seems like Holme's solutions, given in reference 13, 
are the only ones available. These are closed-form incompressible flow solutions for 
rectangular and triangular wings with parabolic arc cross-sections at zero incidence. 

Of the various numerical methods available for the evaluation of ^(x,Y,Z), e.g., 
reierences 14 and 15, the calculative technique employed in the present study is based on 
an application of the method of singularity (ref. 16). Since the superposition principle 
is valid in treating the linear equation (7), this method of solution utilizes a planar 
distribution of elementary sources and vortices to represent the physical wing configuration 
and the induced circulatory system. 
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Thus, for a constant source distribution of strength q(§, C) over an area element AS 

2 1/2 

in the xZ-plane where Z = (1 - M ) z one is faced with an integration of the form 

* CO 

(ref. 16) 


= ' 2 \ f f q(§,C) 2~ 

2TT J J [(x -o 2 + 


w. 2tt / I 2 2 2 1/2 ^ = 

J J [(x -§r + y 2 + (z - c) j ^ 

AS 

Similarly, a constant strength vortex distribution y( 5, C) leads to the representation 


f ^.a 


(x -?) 


y 2 + (z - C) 2 [(x - ?) 2 + Y + (Z - C) 2 ] 


2A72 


d?dC + 


//*• 


0 -o * 2 d?dC 

y 2 + (z - C) 


Since the wing planform geometry is necessarily discretized by a number of small planar 
elements, the appropriate limits of integration are defined by the corner points of the 
individual panels. Within each of these panel elements the strength of the singularities 
q(?/C) and y(?,0 are assumed constant and a control point is selected at which the local 
boundary conditions are to be satisfied. 

The associated computer program which implements the solution of the spherical 

harmonic m = cp + cp of equation (7) for given boundary condition is listed in Appendix B 
s v 

as subroutine CARM . It should be noted that the original listing of this program as obtained 
from the NASA Langley Research Center has been modified to consider perturbation 
velocities in the x-direction only. 


The Integral Equation Formulation 

Equation (2) resembles a Poisson differential equation for which potential theory admits 
the solution 

CO 

cp(x,Y,Z) =cp(x,Y,Z) +cp'(x,Y,Z) - q Q d * dr ' d< * (10) 

3 

— X 

2 2 2 1/2 

and where R^ = C(x - §) + (Y -q) + (Z - Q) ] designates the inverse of the singular 

fundamental solution for the domain of integration. To the harmonic solution cp is added 
another harmonic solution cp'which shall serve as a correction potential to the integral term 

of equation (10). Since the boundary condition (5b) shall also be satisfied by the potential 

cp , the need for the correction term will become apparent when one differentiates 
equation (10) with respect to Y and rearrange to yield 
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00 


< M X,Y,Z > '’V (x ' Y ' Z) =^(x.V,2 ) -±fff C* ? „ 5; ] ? /T1#C ^WdC 

-00 ' ' 

where the kernel 

(Y -tQ 

SY \3 / [(x -§) 2 + (Y - ri) 2 + (Z - C) 2 ] 3/2 

is an antisymmetric function with respect to the xZ-plane for Y s 0, 

. The derivation of equation (10) for symmetrical flows, i .e . withcp'(x, Y,Z) =0, has 
been given by various authors including Gullstrand (ref. 17) who gave a detailed derivation 
in connection with the Green's theorems. Similar derivation is found in reference 18 and 
Klunker (ref. 19) generalizes the use of the integral formulation of the transonic small- 
disturbance velocity potential for the far field boundary condition in connection with 
finite -difference calculation methods, see e .g . reference 20. 


Next, we seek to modify equation (10) for a necessary numerical analysis. If a 
perturbation velocity vector in the transformed space x,Y,Z is defined as (U,V,W) = 
grad cp and if we restrict ourselves to seek a solution for U(x,Y,Z) in the planform 
plane (Y = 0) then a differentiation of equation (10) with respect to x will yield after 
some manipulation 


U(x , 0, Z) - I u 2 (x , 0, Z) = U(x , 0, Z) + cp;(x , 0, Z) + 


■m 


(?,T1,C) 




£-jd§dridC 


O’) 

To facilitate the solution of this nonlinear integral equation, an exponential decay in the 
Y-direction will be assumed for the perturbation velocity component U = U(x, Y,Z), i .e . 


U(x,±Y,Z) - U(x,±0,Z) exp{^Y/r(x,±0, Z)} 


( 12 ) 


where the parameter r (x,Z) = r(x, ±0,Z) is approximately derived from equations (4) and 
(5) as 


r (x, Z) - abs 


(1 -M 2 )'/ 2 

00 

u(x,±0, z) - u 

CO 

/£d 

U 


CO - 

/ dx 


Equation (12) substituted in the integral term of equation (1 1) will then yield, after an 
integration in the r|-direction , the reduced integral equation 
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( 13 ) 


U(x,0,Z) -Iu 2 (x,0,Z) = U(x,0,Z)+cp;(x,0,Z) + ^ 


Jj U 2 (S,+0,C) E + d§d£ + 

ff u 2 k,-0,Q£~ ttdC 

■ — CO 

where the kernel E - E((x-5( ,( Z— Cl ; r ) in the double integrals reads 


E(l x _ ?l / i Z - C|; r± ) “ 2 f e xp(-2Vr ± ) -^"2 (^ _ ] d ' n 

3 § ' 3/. 


n= 0 


?lr / e xp(- 2 Vr ± )|(^-) d n 

n=0 


(14) 


Inserting the relation (valid for Y = 0) 


a n _\ = 


35 W [(x -§) 2 + T ! 2 + (Z-C) 2 ] 3/2 


in equation (14) and making the substitution 


t - r \/\ r i .e. dq = rdt 


where the upper superscripts to the parameter r have been dropped one obtains the 
following expressibns 


E(|*-?I,|Z-C!;r)=|r{ 


x.-5 


a 

/ 


exp(-2t) 


1 


t=0 




372 


dt 


_ 3_ 

3? 




(15) 


Here 0-2 
see e.g. reference 21 . 


M-M 


1/2 


and the integral is obtained from standard tables, 
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rrom reference ll one can aeauce rne following properties 


iH,(a) =H o (a)-lH,(a) 

^N | (a) = -N 2 ( 0 )+jN 1 ( ff ) 

for the Neumann function N and the Struve function H of both integer order v and with the 

v v 

aid of these relations the final evaluation of equation (15) will yield (see figures 1 and 2) 

2 i 

E(|x -§|,|Z - C|;r) “ 4( ^T 1 ) •’ + 

r a La ' ' J L 

t ;( i r) J [ H . w -« 1 ( H i H + N i w ) tN ! ,,) ] ° 6) 

which reduces to 

C7> 

E x = H Q (a) - I (H^a) + N,(a)) + N 2 (o) ] (18) 

for x - §s0 and Z - ( hO respectively. Equation (17) has previously been obtained in 
reference 5, but omits the function r in the denominator as a result of a typographical error. 

Correction for lifting, continuous flows . - In the case of lifting flows, the integral 
terms of equation (13) will give rise to a downwash in the plane of the wing planform. Since 

the potential cp is already required to satisfy the specified boundary conditions in this plane, 

i .e . , cpy = cpy at Y = ±0, the correction potential cp 1 is added to the harmonic solution cp to 

counterbalance this induced downwash. Furthermore , cp 1 must also be a harmonic function, 
but we can confine this requirement (or property) to the xY-plane since equation (13) does 
not depend on any spanwise (i ,e . , approximately the Z-direction) derivatives of the function 
cp' , We will first consider continuous flow and define a correction potential cp 1 =cp'(x,Y,Z) 
which shall be antisymmetric with respect to the planform plane Y = 0 as 

00 

cp'(x,Y,Z) = -I j v'(?,Z) tan" 1 ^—-d? (19) 
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The vortex strength y'( x ,Z) is necessarily defined as 

Y'(x,Z) =i[U(x / +0,Z) - U(x,+0 ( Z)-U(x,-0,Z) + U(x / -0,Z)] (20 ) 

and a differentiation of equation (19) with respect to x together with the limit value 


Y'(x,Z) = l 


im 


Y-*0 it h ' 


Z) 


2 2 

(x -?r + 


d§ 


yields the identitycp '(x,0,Z) = y'(x,Z). S ince a unique circulation is already implemented 
in writing equation ^20) no further considerations have to be made regarding the Kutta 
condition at the trailing edge . 

Correction for lifting, discontinuous flows. - With shock surface discontinuities in the 
flow field a modification of equation ($0) is necessary in order to assure the harmonic nature 
of cp'(x,0,Z) and, hence, enforce a smooth continuous derivative cp'^ . However, instead of 

defining the function cp' we define its gradient in the x-direction as 




(x,0,Z) = ^l%^[U 2 (x,+0,Z) - U 2 (x, -0, Z) ] 


( 21 ) 


and let m(x,Z) represent a spanwise-varying scaling parameter which shall depend on the 
local circulation around the wing in order to satisfy the Kutta condition. Hence, one can 
write 

GO 

m(x,Z) = 2 f [U(',+0,Z)-U(?,+0,Z)-U(?,-0,Z)+U(?,-0,Z)][U 2 (§,+0,Z)-U 2 (§,-0,Z)] _1 d§ 


( 22 ) 

It should be emphasized that equation (21) is not unique and that a better downwash repre- 
sentation at the leading 
factor [1 + (9y /dx)^]' 


sentation at the leading edge could be achieved with, say, an adoption of the Riegel's 


From the elementary relation 

U 2 (x,+0,Z) - U 2 (x,-0,Z) = [U(x,+0,Z) - U(x,-0,Z)][U(x,+0,Z) + U(x,-0,Z)J 

where the first and second factor on the right hand side represents the linearized perturba- 
tion velocity due to lift and thickness respectively. One sees from equation (21) that the 
distribution in chordwise direction of the x-component of the correction potential is defined 
to be of a similar form as the disturbance due to thickness. This distribution, however, is 
modified by an amount proportional to the local chordwise lift. It should be noted that the 
use of a different combination of the linearized velocity components in equation (21) could 
be better defined from known lifting flow calculations. 
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NUMERICAL ANALYSIS 


The Algebraic System 

The solution to the nonlinear integral equation (13) cannot be obtained analytically 
for wing configuration of arbitrary planform and thickness distribution. Hence, a numerical 
treatment of the general problem is mandatory and the first step in this direction is to 
discretize the area of integration in accord with the division of the wing planform as used 
for the numerical evaluation of the harmonic' solution . The range of integration shall 
be confined to the wing planform and the nonlinear compressibility sources outside the 
cylinder which encloses the planform and which have the generatrix parallel to the 
Y-direction are therefore not being considered. 

Let the upper and lower surface of the left (or right) half of the wing planform each 
be represented by N number of trapezoidal panels and let the unknown function U = U(x,±0,Z) 
assume a constant value within these specified panel boundaries, see figure 3. Then equation 
(13) can be replaced by a system of 2N algebraic equations which are nonlinear (i.e., 
quadratic) with respect to U = U (x,0,Z), i.e. 


2N 


1 . .2 rr 


F i< U l' U 2 U 2N> = U i "2 U i - U i +C i + Z e ii U i =0 


i=i 



. ,2N (23) 


where the index i ^ N shall identify a value on the upper surface of the wing (or panel). 
Consequently, the lower surface values are associated with the range N< i< 2N. 

Due to the lateral symmetry assumed for the wing planform and its position relative to 
the oncoming flow (i .e . wj= 0) only one half of the wing need to be considered in the 
algebraic system (23) since the influence from two panels which are positioned symmetrical 

with respect to the x-axis can be lumped together as e.. = e.. + e.. . (Note that U. = 11. 

•I -I, <l 2 ' 2 

where the indices 1 and 2 refers to each half of the wing.) 


The influence coefficients e . j a re integral representations of the influence function 
E = E ()x - ?| , | Z - £|; r) as evaluated over a surface element AS , i.e. of the form 


4rr 


Ih 

AS 


?|,|Z -C|;r)d?dC-e.. 


for | x — ?| ¥■ 0 and |Z 0. Due to the particular and perhaps unfortunate influence 
pattern of positive and negative contributions from the integral terms of equation (13) 
care must be exercised in the numerical evaluation of the coefficients e.., see figure 4. 
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Th is is further illustrated by the series representations 


E (p) = -4 


E Z ( X ) - 4 


1 2 4 2 

~2 "" 3 P + 45 0 


1.1 1 3. 

X 3 X_ 45 X + 


|p|-0 

(24) 

lxl-0 

(25) 


of the influence function E in chordwise and approximate spanwise directions respectively 
where p = 2j Z - Cj/r and x = 2jx - ?|/r . Since both equations (24) and (25) show 
similar behavior at the singular points p= \ = 0, but of opposite sign, we let e..= 0. 

Furthermore, we define the influence coefficients e.. as 

'I 


e .. = ^ [E (') + E (2) + E (3) + E (4)3 


(26) 


Equation (26) is derived from a simple application of the mean value theorem of integral 
calculus by relating the influence strength e.. at the centeroid of each panel to the 

influence measure at four points 1= 1, 2,3,4 located on the four panel boundary lines 
(see figure 5). A better approximation of e.. would be achieved, however, if a surface 

11 (4) 

spline function is utilized in connection with the four values E' and the local value of E at 
the centeroid, see e.g. reference 23. Note that for symmetrical flow problems it is sufficient 
to consider only a system of N equations with the minor modification of doubling the value of 
the influence coefficients e..(i, j = 1 ,2, ... ., N). However, the following analysis will be 

presented for a system of 2N equations. 

The function C. (i = 1 ,2, . . .,2N) appearing in equations (23) is a velocity correction 
function associated with circulatory flows and is easily obtained from equation (20) as 


C. = A A. 

i 2 i 


i = 1,2, ...,2N 


(27) 


for continuous flow. For circulatory flows with a shock discontinuity, equation (21) approxi- 
mately defines the required correction as 


N 


1 


C. = =■ B. y (1 -6..)=i 
I 2 I II B. 


i =1 


i = 1,2, ...,2N 


(28) 


where 6.. represents the Kronecker delta and where 
M 
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u. - u. - u.,. . + U., KI 
i i i+N i+N 

i ^ N 

u. - u. - u. . . + u. Kl 

i i i-N i-N 

i > N 

-2 —2 


U i ' U i+ N 

i < N 

-2 —2 


ur - u: . . 

i i-N 

i > N 


In the case of non-circulatory flows, equations (27) and (28) will yield the expected 
corrections of C. = 0. System (23) is now completely determined and the task of finding 


a vector U.(i = 1,2,... ,2N) which satisfies equations (23) for given £.. and U. wil 
the subject of the following section. 


be 


Methods of Solution 

Solutions to the nonlinear system (23) could be obtained by various iteration or 
relaxation schemes and each method has its own merits and limitation. The decisive factor 
in selecting a solution technique, however, is the type of flow under consideration, see 
figure 6. Since we are also interested in obtaining discontinuous solutions which shall 
represent flows with compression shocks, we will distinguish between three different 
types of flow problems, namely 

o Subcritical flow 

o Supercritical flow without discontinuities 

o Supercritical flow with discontinuities 

and describe some numerical approaches for their unique solution. In essence, each one of 
the numerical methods to be introduced should yield the same solution to the quadratic 
equations (23) for identical given flow configurations. This is especially correct in the 
subcritical case where only the smallest roots of system (23) are sought (i .e . , all U. -values 

corresponds to subsonic ve locities) and this single solution is necessari ly unique . 

However, a supercritical discontinuous solution vector U. (i = 1 ,2, . . . ,2N) of 
system (23) is signified by being a combination of roots which does not belong to the 
same family, i ,e . subcritical or supercritical where the critical value U.= 1 correspond to 

the case of two equal roots . Hence, multivalued solutions are possible and some additional 
condition must be introduced in the calculation scheme in order to render the obtained 
solution unique. Such a condition is in our case a substitute for the entropy inequality 
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which forbids the occurrence of expansion shocks in real flows, and the present study has 
adopted a version of the method of steepest descent to be a suitable calculation scheme 
in which the "entropy" condition can be introduced. This is achieved by an approximate 
predetermination of the location of the critical points (U. = 1) as indicated by the linearized 

solution U. (i = 1 ,2, . . .,2N) and this again forces the calculation scheme to search for 

one particular solution. 

An intermediate class of transonic flow solutions are characterized by being both 
continuous and supercritical in the sense that subsonic and supersonic velocities can 
exist together within a certain control volume without being separated by a surface of 
discontinuity. Both theoretical and experimental investigations made in the more recent 
years have confirmed this type of flow phenomena and the introduction of the method of 
parametric differentiation in the present study can be viewed upon as an attempt to 
calculate and to plausible explain such flows. This has been achieved by embedding the 
solution to the system (23) in a Cauchy type problem in which a value of U. = 3/2 appears 

as the upper critical value for a unique dependence of the initial data upon the sought 
solution . 


In summary, it should be emphasized that the economics of finding a numerical 
solution of the stated system of equations (23) has been considered in the selection of a 
suitable alogrithm. Thus, the fast converging method of Newton is found to be far the 
best choice for a subcritical flow calculation despite the fact that the other two solution 
techniques could be applied as well. These two methods, on the other hand, represent 
choices for the calculation of the more intricate supercritical flow cases and the ambitious 
user will find ample opportunities to make improvements of the listed computer program. 


The Newton's method. - An application of Newton's method of successive approxi- 


mation to the solution of system (23) will generate a sequence of vectors U. 
the relation 


£k+i ] 


from 


u( k+ I} _ yW - [F^OJ .)]" 1 F^(U.) 


= 1,2, ...,2N 


(29) 


where the integer in the bracket { } identifies the sequence number and where the 


starting vector U. 


{k=0} 


can be chosen to be identical to the known vector U. . This 

,{k] 


method is particularly valuable for subcritical flow calculations , i .e . U ! J < 1 

(i = I , 2 , . . . , 2N), for which the guarantee of uniqueness of a solution be longs to a 

classical proof in fluid mechanics. (See e.g. reference 24 and its list of references). 

As long as the functional determinant (or Jacobian) is nonvanishing and positive for 

al I values U. < 1 , the necessary inversion of the Jacobian matrix J(U.) = F. , (U.) or 
i i U. i 

i 

J(U 1 ,U 2 ,...,U 2 n) = 5( F 1 ^ F 2'-’ ' '^2N^ / '^2' ’ ’ ‘ '^2N^ 

can be handled by numerous available computer subroutines. 


15 



Method of parametric differentiation . - A continuous application of Newton's 
method to the solution of system (23) for values of U. exceeding unity (i .e . for super- 
critical flow), would imply some arbitrariness with regard to the uniqueness or the 
usefulness (expansion shocks can occur) of the obtained solution. For our purpose in 
controlling the development of an iterative unique and continuous solution to the system 
for the case of U.> 1 we will adopt the method of parametric differentiation (see e.g. 

reference 25) which has found wide applications in mathematical analysis, including 
a treatment of the transonic flow problem in two dimensions (ref. 26). 

Our method of solution will be based on the particular one-parameter operator 
homotopy (embedding) 


G.(n;U r U 2 , 


...,U 2N )=U. -U. + p[F.(U) -U. + U.] 


= 1,2, ,2N (30) 


where pi is a real-valued parameter in the range 0^p£ 1 . From the implicit function 
theorem as applied to equations (30) one obtains the differential properties 


dU. 

i 

dTT 


= -J"V;U 1; U 


2 ' 


BG . 


,UokJ 


2N Bp 


i = 1 , 2, . . . , 2N 


(31) 


where the Jacobian matrix of the system is 


J(l_i; Uj ,U 2 , . . .,U 2N ) B(Gj ,G 2 , . . .,G 2 ^) / d(U^ , U 2 , . . . ,11^) 


The first-order nonlinear ordinary differential equations (31) together with the 
initial conditions 

u = 0: U. = U. i = 1,2 ,2N (32) 

defines a Cauchy problem which at p = 1 yields our sought solution of the algebraic system 
(23). There exists various powerful methods for the numerical solution of system (31), see 
for example reference 27, and we will confine ourselves to the use of Hamming's modified 
predictor-corrector method together with a fourth order Runge-Kutta method for the compu- 
tation of the starting values. Since the right-hand side of equations (31) includes the 
parameter p, the system (31) was solved at every appropriate stage of the integration 
(0 ^ p <■ 1 ) by including a modified Gaussian elimination procedure . 


A condition for a unique solution of system (31) in the range 0 s p s 1 is the non- 
vanishing of the functional determinant (or Jacobian) of G. (p; U^, . . ., U^) w '^ 

respect to U. for fixed p . This determinant is positive under the condition that the 
elements on the principal diagonal are positive, i ,e . from equations (23), (27), and 
(30) we obtain 

1 - p( u i - ^) > ° i = 1 ,2, . . . , 2N (33) 
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since rne initial solution U. Is necessarily continuous and will yield a positive 

Jacobian, one can conclude by virtue of inequality (33) that a unique and continuous 
solution of system (31 ) can at least exist for U. < 3/2 (i = 1,2,..., 2N) . 

However, it is more consistent with the present analysis to express this limit for 
continuous flow directly in terms of a maximal local perturbation velocity 



3(1 -M 2 ) 

CO 

2(x + 1)M 2 

CO 


and then apply, for example, the isentropic relation M = M(x,u/u , M ) for the evaluation 

00 CO 

of the corresponding Mach number. The dependence between this maximal local Mach 
number and the AAach number at infinity is given in figure 7. It is seen that the predicted 
limitation of shock-free transonic flow is in general accordance with the two-dimensional 
results from Nieuwland's theory (ref. 28) and of Korn (ref. 29) depsite our underlying 
assumption of small disturbances. Also depicted in the figure is the stability criterion of 
Spee (ref. 30) which limits the value of the local Mach number to (5/2) ' forx = 7/5. 

The method of steepest descent. - Fora certain value of the parameter p the 
Jacobian matrix in system (31) might become singular and the solution curve to the 
Cauchy problem, equations (31) and (32), can split into two or more solutions . Such a 

point p.= p* in the interva I 0 < p -s 1 for which det J(p.; , U^, . . . , ) = 0 can be 

regarded as a bifurcation point and the associated solution Ut = U.(jj*)will be designated 

the bifurcation solution. Since the next step in the integration of system (31) would 
involve the crossing of a line of a singular Jacobian matrix, a continuous dependence of 
the solution on the initial data is no longer possible. The condition that the bifurcation 
solution U*satisfies equation (30) could be utilized in connection with the introduction of 

a family of differentiation parameters in order to continue the solution curve past the 
singular lines. This difficult problem, however, will not be further pursued in the present 
report. Instead, we will apply the method of steepest descent to the solution of system (23) 
which represents a set of quadratic equations with respect to U.(i = 1,2,.. . ,2N) and the 
roots of the system for the case of discontinuous flow are easily obtained as 


2N 


1/2 


U. , = 1 + sgn 

' 1,2 


1 - 2 


U. - C. - 2 


Y e.. uf 

'I I 




i = 1 , 2, . . . , 2N (34) 


where sgn - ±1 depends upon U.^; 1 . Hence, a bifurcation of the solution to equations (34) 
subject to the initial data U. = U. occurs when the discriminants vanish and this is equivalent 



to considering the identity of the relations 

2N 




i = 1,2, ...,2N 


(35) 


Furthermore, the inequality of (35) gives the condition for positive discriminants and hence 
the condition for real-valued roots of the algebraic system. 


Let a function iji - f(f^ defined as 

2N 


^ S V U r U 2' * " ,U 2N^ 


wh 


ere 


f.(U, ,U 0 , . . -/U 2 n) = 1 - U. + sgn 


1 


u. -c. - 

t I 


(36) 


2N 

2 S e .. U 

^ i| | 


i =1 


1/2 


i = 1,2, 2N 


(37) 


which is obtained from a rewriting of equations (34). In order to find the roots of system 
(37) the method of steepest descent considers the variational problem of searching for the 
zero value of i|r from equation (36). The calculation scheme in the method produces a 

sequence of vectors [U. } from the relations 


ai(f 00 fW f tk}j 
ylk+l] = u w _ {k} an 1 ' 2 -' f 2N } 

i i BU. 


i = 1,2, ...,2N 


(38) 


where the step length \ along the gradient direction of the function i|f is determined at each 
sequence number k as 


00 - 


, 00 - 


2N r 

X 




{k} 


dU. 


The sign convention in system (37) is either determined from the condition 

(k] 

U. $ 1 : sgn - ± 1 or from some information which yields the approximate geometric 
location of the sonic line and the foot of the shock discontinuity. 
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The present method utilizes the linearized solution U. (i - 1 ,2, . . ,,2N) to 

approximately define any supersonic region on the wing surface or in other words to 
enforce the sign convention in equations (37). This is accomplished by mapping the 
region for 

(r+ 1)M 2 

- — CO 

U. > yj- C* 

' 1 - M 2 P 

CO 

* 

where Cp is the critical pressure coefficient as obtained from the exact isentropic 

relation and then shift the upstream boundary line one chordwise panel length downstream. 
The resulting region is then our definition for the area of supersonic velocities and a 
positive sgn-value in system (37) shall only be connected to points U. which lies within 

this region. An improvement of this rather arbitrary technique of defining a sonic line 
and a line of discontinuity has been attempted with a partial use of the associated 
bifurcation solution , but required computation times have made such an improvement 
impractical . 

RESULTS AND DISCUSSION 
Two-Dimensional Flows 

As first examples of continuous flow, subcritical pressure distributions past a 
nonlifting and lifting NACA 0012 airfoil have been calculated using the Newton's method 
and are shown in figures 8 and 9, respectively. Comparisons are made with "exact" 
numerical solutions, see reference 31, which are obtained from an application of Sells 
method to the full two-dimensional gasdynamic equation. Furthermore, it is found that 
the present transonic small-perturbation solution as seen in figure 9 is very similar to the 
solution given by Krupp (ref. 32) which also assumes small disturbances to the freestream 
condition. 

Figures 10 and 11 show results for supercritical flow past the same airfoil configurations 
as given in the two preceding figures, but the present solutions are continuous in the sense 
that supersonic velocities exist on a portion of the profiles without being terminated by a 
shock. These results are obtained by the method of parametric differentiation and are 
compared with solutions which have been calculated with the Jameson program, see e.g. 
reference 33. It is seen that these latter solutions exhibits small compression shocks, at 
least in the solution given in figure 11 . For freestream Mach numbers slightly higher than 
those given in the figures, no continuous solution could be obtained. This occurrence is 
interpreted as the Mach number limit for which shock discontinuities start forming on the 
airfoils . 

The result for discontinuous flow past an 8.4 percent (t = 0.084) thick symmetrical 
parabolic arc airfoil at zero angle of attack is given in figure 12 and compared with the 
numerical solutions of references 20 and 34. The solution associated with the latter 
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reference is, for the given thickness parameter, related to a freestream Mach number 
of = 0.8496. The lifting case of a discontinuous solution is depicted in figure 13 

for a NACA 0012 airfoil and comparison is made with the numerical results of reference 
35. 


Three-Dimensional Flows 

To test the Woodward/Carmichae! computer program which calculates the linearized 
solution to the posed wing problem, a subsonic lifting case for a rectangular wing of 
parabolic arc cross-section has been evaluated and plotted in figure 14. The chordwise 
pressure distribution at various spanwise stations (identified by the spanwise coordinate C 
which is normalized to the semispan) are in good agreement with the results of figure 10 
from reference 7'. The familiar singular behavior of the solution at the leading edge 
region, however, is apparent. This is also brought out in figure 15 which compares the 
calculated incompressible flow solution at mid-span for a high aspect ratio wing with a 
blunt leading edge with the exact two-dimensional data (reference 36). A further test of 
the accuracy of the linearized solution is depicted in figure 16 which compares the 
numerically obtained results for a symmetrical wing with the analytical solution of 
Holme (ref. 13). It is somewhat surprising to see that the numerical solution yields 
higher perturbation velocities than the analytical. The reason for this might be in the 
low number of panels (i .e . , 100) used to represent the semi-wing planform in the 
numerical calculation scheme. 

Figures 17 and 18 compares the present subfcritical results, as obtained from an 
application of Newton's method, for a high aspect ratio non-lifting and lifting wing 
respectively with the "exact" two-dimensional values given in reference 31 . For this 
aspect ratio of 7 the flow at the midspan of the wing should be nearly planar and figure 18 
also partly validates the downwash correction potential (see equation 19) employed for 
this case. The sensitivity, however, of the number of chordwise and spanwise panels used 
in the calculation scheme is brought out in both figures. Similar results for the same 
non-lifting wing, but of smaller aspect ratio, is given in figure 19. 

For an arbitrary wing configuration the planform of a trapezoidal half-wing is 
defined by its four corner points and the profile geometry is input in a table-form. This 
is illustrated by the swept wing configuration of figures 20 and 21 . The cross-section of 
this wing is a RAE 101 airfoil and the profile surface coordinate, slope, and curvature at 
various chordwise stations are tabulated in subroutine SECTIN as ZTAB1, ZPTAB1, 

CUTAB1 , and XTAB respectively. Calculated supercritical results were obtained by the 
method of parametric differentiation and compared with, the experimental data of reference 
35. Both types of data indicates that the flow solutions are continuous. 

Discontinuous flows are obtained with an adoption of the method of steepest descent 
and sample flow calculations are given in figures 22 through 25. The first two figures 
show results for a non-lifting parabolic arc wing at two different spanwise stations. 
Comparisons with the numerical solutions of reference 7 are not very encouraging and this 
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is even more true for the lifting case as shown in figure 24. However, if the dis- 
continuous supersonic region is confined to only a small part of the wing (e .g . in the 
root section area) a better agreement of the results are found (figure 25). 

CONCLUDING REMARKS 

The 1 compressibility effects on the flow past lifting and nonlifting wings at high 
subsonic Mach numbers have been studied with an application of the transonic integral 
equation method. This semi-analytical calculation method preassumes a solution for 
the equation describing linearized flow around the same wing configuration, i .e . assumes 
the availability of the PrandtUGIauert solution, and this problem is solved by means of 
known analytical or numerical solutions. The aim of the study was to incorporate the 
nonlinear compressibility effects in a pressure calculation scheme which also should include 
the case of a shock compression surface in the flow field. 

The present calculation method involves a summation of some compressibility effects 
which are distributed over the wing planform area and which are a function of local 
velocity and profile curvature. Since these nonlinear effects can attain different sign 
convention an approximate numerical integration over the domain of influence depends 
largely on the grid size employed. An alternative integration technique (approximate 
two-directional) which involves one sign convention for the influence measure only, 
has therefore been developed and applied to all calculation cases shown in this report 
(see figure 3). 

The practical desirable case of shock-free supercritical flow was studied by an 
incorporation of the method of parametric differentiation in the numerical evaluation of the 
governing algebraic system. Obtained results show that local supersonic velocities can 
indeed exist in the flow field without a discontinuous recompression in the form of a 
shock surface. It is felt that these findings constitute an important contribution to the 
subject involved . 

For the calculation of supercritical flows with shock discontinuities, however, the 
present methods seems less effective. The reason for this lies in the fact that points on 
the sonic line represent extreme values for the stream-density quantity and an accurate 
influence definition in the neighborhood of these points is required. Also, a different 
iteration scheme for the evaluation of the nonlinear equations could greatly improve the 
usefulness of the presented analysis technique. ; 
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APPENDIX A - BASIC EQUATIONS FROM TWO-DIMENSIONAL ANALYSIS 


For wings of infinite aspect ratio the analysis becomes two-dimensional and available 
results for this case can be utilized in the calculation of flows about finite wings (see 
reference 5). The appropriate integral equation (as compared to equation (10) on page 7) 
will then read 


cp(x,Y) - cp(x,Y) + cp'(x, Y) - 2^- Jf [cp ? 


^sW n iq d5dT1 


(Al) 


2 2 1/2 

where R 2 = [(x - §) + (V - rq) ] ' . A similar assumption for the decay of the perturbation 

velocities in the Y-direction as made in the foregoing three-dimensional analysis will reduce 
equation (Al) to an equation involving the unknown quantity U = U(x,±Y) of the form 


CO 

U(x,0) - 1 U 2 (x,0) = U(x,0) +cp^(x,0) + f U 2 (§,±0) E 2 (|x - § | ;r ± ) d? (A2) 

— 00 


where the two-dimensional influence function - E 2 (]x -§|;r) is given as 


E 2«- E 2 i^) 

sin (x) - Si(x)J - cos(x) Ci(x) 


4 

tt 


(A3) 


Here Si and Ci are the sine and cosine integrals respectively. 

An analogous numerical integration of equation (A2) would involve the determination 

of influence coefficients e.. and this leads to a definition of 

'I 

Ax/2 

e ii = 2F / e 2 d? 

0 
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which yields after a substitution of equation (A3) in the integrand of equation (A2) the 
following value 


E ii ( \> ) = E ii(r) 

= i- cos(x o ) [ Si(x o ) - y] - S in(x o ) C i<x o > + J 

The values of the coefficients g j j (where j / i), however, is determined direct from the 
influence function E 2 as e I j = AxE2(|x - 5|;r)/(4r). These influence coefficients are used 
in the approximate integration of the nonlinear influence measure over the wing, see e.g. 
figure 3 . 
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APPENDIX B - LISTING OF COMPUTER PROGRAM 


A computer program has been written in the FORTRAN IV language to implement the 
numerical calculation schemes described in the foregoing sections. The program consists of 
a main calling program designated MAIN and 17 external subroutine or function subprograms. 
Subprograms which are part of the particular Carmichael computer program received from the 
NASA Langley Research Center are designated by an asterisk (*). A flow chart showing the 
calling sequence of the program subroutines are given in figure 26. 


Name 

MAIN 

CYL 

SICI 

ACSH 

GAUSS 

HPCG 

FCT 

OUTP 

CARM* 

WNGEOM* 


Description 


Main program 

A subroutine subprogram which evaluates the cylindrical functions N 
(Neumann function or Bessel function of the second kind) and H ' 
(Struve function) both of integer order v. 

A subroutine which evaluates the sine and cosine integrals (Si and Ci) 

A function subprogram to evaluate the transcendental function y = arcsinh x 
(arcus sine hyperbolicus). 

A subroutine which solves a system of linear algebraic equations by a 
modified Gaussian elimination method. 

A subroutine to solve a system of first-order ordinary differential equations 
A subroutine which defines the differential system (31) 

A subroutine used for control of output from HPCG 

A subroutine used as the calling program for the evaluation of the arbitrary 
harmonic or linearized solution (i.e. a modified Carmichael program). 

A subroutine to evaluate the geometric characteristic of the wing. 


FILL* 


A subroutine for function interpolation in connection with WNGEOM. 
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SECTIN* 

TAINT* 

EVAL* 

COMP* 

(TCOMP) 

INVERT* 

FORCE* 

VECON 


A subroutine used to input the profile geometry. Note: The four data 
tables XTAB, ZTAB1, ZPTABI, and CUTABl given in this subroutine are 
for a RAE 101 airfoil and they represent 19 chordwise values for the profile 
chordwise station together with the associated measure for the thickness, 
slope and curvature with reference to a profile thickness ratio of one. For 
arbitrary profile geometry appropriate input data tables must be supplied by 
the user. If the number of chordwise input stations are different from 19, 
a dimension statement and three sequential calling statement for subroutine 
TAINT must also be changed accordingly. 

A subroutine used in connection with SECTIN to interpolate geometric data. 

A subroutine to evaluate the aerodynamic influence coefficients in connection 
with the harmonic solution. 

A subroutine used for the calculation of the linearized perturbation velocities 
due to thickness. This subroutine has an alternative entry point identified by 
the name TCOMP. 

A subroutine for the matrix inversion of the influence matrix evaluated in 
EVAL. 

A subroutine used for the calculation of the complete linearized solution in 
terms of perturbation velocities. 

A subroutine which converts local perturbation velocities to local Mach 
numbers, pressure coefficient and pressure ratio. 
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Listing of main calling program and associated external subroutines. 


c main program 

c 

c this program calculates the transonic pressure distribution 
C ON arbitrary lifting WINGS at subsonic speeds 
C Based ON A (NONLINEAR) integral equation method. 
c 

C The PROGRAM was LARGELY DEVELOPED UNDER CONTRACT NASl-10665 

c sponsored by the loads division 'of the nasa langley , 

C RESEARCH CENTER. BUT INCORPORATES ThE WOODWARD/CARMICHAEL 

c Computer program for the calculation of the linearized solution. 
c 

c the program is dimensioned for a total of loo panels 
c To represent the wing PLAnFORM and this LIMITS the number 

C OF PANELS WHICH CAN Be USED FOR ThE LIFTING CASE TO 50. 

C A SYMMETRICAL NONLIFTING WING. HOWEVER. CAN BE REPRESENTED 

C BY A MAXIMUM OF 100 PANELS. THE STORAGE REQUIREMENT 

C IS APPROXIMATELY 145 COU IN OCTAL. 

C 

DIMENSION 

1 .lUX(lOC). IUZ(IUO), UllOO). U0LC100). FINT(IOO). 

2 DFU(IUO). PRMT ( 5 ) » EP(4). DELTX ( 50 ) . AUX ( 16 . 100 ) ► 

3 XSU ( 50 ) * XSL ( 50 ) « XXU<50>» XXL(50) 

C 

common /dy/ dery(ioo) .ycioo) »ui.(ioo) *dp 
Common /hol/ eps<iooOo)» a(ioioo) 

Common /prndtl/ upgUgo) .ilift.nrun.symf ; . " 

common /vel/ mach, mac,, su.betasq. beta. aami.aam2.aam3 
Common /para; is/ nwing , panels . sref , refmom , cbar . span . oc 
Common /wngpms/ rootu> .tiP(4) .m.n.type.F(ioi) ,G(iod .p(ioi) . 

1 SHEAR ( 101) , ISECT.THIcK(b) 

COMMON/XE.V/XW ( luO ) . ZW ( 100 ) 

Common /cukv/ cuRdooj 
REAL MACH, MACHSQ.BETASQ. BETA 

namelist /flo*/ hach.gelm.ndelm 

NAMELIST /wing/ iwing, influ, method .root. tip, thick, ISECT. Lx. lz» 

1 SPA,,, T AO 

namelist /shock/. xsu.aSl.xxu, xxl 

C ' • ' 

c constants and the input 
c 

IT = 1 

ACC=0.001 
Pl-3«l4l5926B.55 
PID=2.0*PI 
P I 4=h . 0 *P I 
Pld3-5.0*PI/3«0 
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non o o o 


PlH=l. 5707963 
Pl 22 = 0. 3163099 
Pl 2 P=G. 5 /PI 
Pl 2 =u « 63 obl 977 
Pl 3 = 3 . 0 *PI 
Pll=i. 0 /Pl 
PI 4 I= 1 . 0 /PI 4 
TyP£=4 
0 C =2 
C 

10 READl 5 »rtiNG) 

IF (IftING.EQ. 9 ) GO TO 9999 

Read ( 5 *flo«) 

lLlFr=0 

IF (IwING.GT.O) go TO 100 
SPAN= 2 . 0 *TIP( 3 ) 

SEMI=SPAN/2.0 

RGT=K00T(2)-R00T(1)+TiP(2)-TIP(D 
Sr£F=TIP(3>+ROT 
ReFmOM= 0 . 26 * (ROOT ( 2 )— ROOT ( 1 ) ) 
CbAR=ROT/2.0 
TaD=TNICM1) 
ar~span*span/sref 
Go TU 110 
100 TaUP=TAU/PI 
S=bPMN/2.0 
SEMI=S 

T A*-‘S = T AU*StMI 
110 MaC,HS 0 =MACH*MACH 
SETASuxI.-MACHSn 
0 ETA=SuRT (ABS(BETASQ) ) 
AaMi= 10 . 0 /( 7 . 0 *MACHSQ) 
AAM 2 = 0 . 2 *MACHSQ 
AAM 3 = 0 . 7 *MACHSQ 
am=mach 
N?LX 
M=LZ 

LXZ=UX*LZ 

LL=LXZ+LXZ 

Ll 2 =LL*LL 

LX 22 = 1 -XZ*LXZ 

0 ELTZ=SEMI/FL 0 AT(LZ) 

OLZ 2 =uELTZ/ 2.0 

nrun=u 


The program header 


200 WRITE (5*210) 

210 Format uhi*27x*36h Velocity and pressure distribution * 

1 40 HON AN ARBITRARY WING AT SUBSONIC SPEEDS /42X* 

2 4faH BASED ON A TRANSONIC INTEGRAL EQUATION METHOD ////) 
WRITE (6*w'lNG) 

wR I f E ( 6 * FLO'.v ) 

The prandtl-glauErt solution 

215 AR2=AM*AM 

UET 2- 1 • 0-AH2 
BETA-SORT IHET2) 





IWl=IWlNG+l 

GO TO ( 460 » 217 » 420 ) » IWi 
217 Ar=2.0*SLMI 

DlTx=1.0/FLOAT(LX) 

Do 200 NZ=1#LZ 
□ELTX(UZ)=OLTX 
Dx2=l)LTX/ 2.0 
Do 220 NX=1»LX 
' K=(NZ-1 )*lX+NX 
zw ( K ) =-DL22+FL0AT ( NZ ) *DELTZ 
X * ( K ) =-OA2+FLOA T < NX ) *oLTX 
220 Continue 
230 Continue 

DO 300 NZ=1»LZ 

Nn=(NZ-1)*LX+1 

SMZ=(S-ZMNN) )*c3ETA 

SpZ=(S+Z*(NN) WuETA 

DO 300 NX=1»LX 

K=(NZ-1)*LX+NX 

X=2.U*XW(NX)-1.0 

XP1=1.0+X 

XMl=1.0-X 

Bl=ACSH(XPl/SPZ) 

B2=ACSH(XM1/SPZ) 

B3=AC5H(XP1/5MZ) 

B4=ACSH(XM1/SM2) 

B5=ACSH(bPZ/xMlJ 

B6=ACSH(SMZ/XM1) 

B7-ACSH(SPZ/XP1) 

B6=ACSHCbMZ/XPl) 

DuM=SPZ*(Bl+t32)+SMZ*(ti3+B4>-X*(B5+B6-B7-B8> 

UpG ( K ) =T AUP+UUM/BET A 

CUF<{K)=-4.0*TAU 

IDX(K)=NX 

IDZ(K)=NZ 

300 continue 
G o TO 480 
420 Ak=4.0*SEMI 

DO 330 NZ=1#LZ 

DlTX=(S-FL0AT(NZ-1)*DELTZ-DLZ2)/S 

DlLTX(n'Z)=DLTX/FLOAT(LX) 

DLX2=DELTx(NZ)/2.0 
do 320 NX=1#LX 
K=(NZ-1)*LX+NX 
Zw(K)=-DLZ2+FL0mT(NZ)*DELTZ 
Xw ( K ) =1 . O-DLTX+FLOAT ( uX ) *0ELTX C NZ ) -DLX2 
320 CONTINUE 
330 CONTINUE 
SsS*BET A 
SS=S*S 
S0 = 1 *0+55 

Sis(1.0+2.0*SS)/(S*S0) 

S2=2.0/S 

S3=2.0*5/S0 

S4=l.O/SGRT(SO) 

S5=l«0/5 
X0F=(1.0-3EX2)/S 
DO 4i>0 NZ = 1»LZ 
NN=(NZ-1)*LX+1 
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Z=ZW(NN)*BETA 

zz=z*z 

SPZ=S+Z 

SMZ=S-Z 

Zl=S.-iZ*SMZ 

Z3=spz*spz 

A6=2.0*Z 
DO 4b0 NX-1 > LX 
K=(NZ-1)*LX+NX 
X=XW(K) 

XX=X*X 

Xl=(1.0-X)**2 

XS2=2.0*S*X 

XZ=SuRT(XX+ZZ> 

Al=(b*X-Z)/S0 

A2=SPZ-XS2 

A4=(S*X+Z)/S0 

A5=S.“IZ-XS2 

SxZM=S*X-Z 

SxZP^S *X+Z 

Vl=(1.0-X + S*Si' / lZ)/SXZM 

V2=(X+S*Z)/5XZM 

V3=(1.0”X+S*SPZ)/SXZP 

V4=(X-S*Z)/SXZP 

Xl 0 = 1 . 0 -x 

V5=SMZ/X10 

V6=SPZ/X10 . , * • 

V7=Z/X10 

B6“ACSH ( Vl ) 

67=ACSH(V2> 

b9=ACSH(V3) 

6lO-ACSH( V4 ) 

Sh5=ACSH(V5) 

Sh6=AC5H ( V6) 

SH7=ACSH( V7) 

B2=SORT(Xl+Zl)+SQRT(Xl+Z3) 

B3=S2*S3RT(Xi+ZZ) 

B4=S3*XZ 

B5=S4*(A1-A2) ' 

88=S4* ( A4-A5) 

Bll=A2*SH5 - • ’ 

Bl2=Ab*SH6 

Bl3=«d*SH7 1 

DuM=bl*a2-b3-B4+B5*<B6+B7)+B8*(B9+Bl0)+S5*(Bll+B12+B13> 

UpG ( n ) =TaUP*OUM/BETA 
CuR(X)=-4.0*TAU 
IDX(K) =NX 
Il)2(K)=NZ 

450 continue 

go TO 430 

460 call carm 

do 470 NZ = 1 f LZ I . 

NW=(NZ-1) *LX+1 

i-JNl=NN+l 

OELTX < NZ ) =Xw ( NN1 ) -XW ( |.N) 

00 470 NXzl » LX 
K=<NZ-1)*LX+I4X 
iOX(i\)::NX 
IDZ(K)=NZ 
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470 continue 

480 CONTINUE 
NS=LL 

IF (IL1FT.EQ.0) NS=LXZ 
NS2=NS*Nb 

IF (ILIFT.EO.0) GO TO 500 
00 490 I=1*LXZ 
II=I+LxZ 
Xic, 1 ( 1 1 ) =Xw ( I ) 

Zw(ID=ZU(I) 

CUR(II)=CUK(I) 

IDX ( I 1 ) =IDX ( I ) 

IDZ ( I I ) =IDZ ( I ) 

IF (1WING.EQ.0) GO TO 490 
UpG (II) =UPG ( I ) 

490 continue 

the reduced variables 

500 NIT=0 

AM7=U.7*AM2 

AM22=0.2*AM2 

AMF::i0.0/(7.0*A;-l2) 

XM=AM*S3RT( 1.0/ <1.0+0. 16666666* (AM2-1.0) ) > 

CpST— 2.0/ (1 ,4*AM2) *( ( (2.0+0.4*AM2)/2.4)**3.5-l.0) 

FaK=2.4*aM2 

FAT=oET2 

ry=beta 

ru=fak/fat 

RV=RU/RY 

UCR=-RU*CPST/2.0 

AMT= ( 2 . 4*AM2 ) **U .33333333/TAU**0 .66666666 
EMR= ( AM2-1 • 0 ) / ( 2 • 4*TAu*AM2 > **0 . 66666666 
HMT=rAK**0 • 33333333 /TaU**0 *66666666 

hmt=hmt/amt 
DO 8U0 1 = 1 * NS 
UL<I>=RU*UPGU) 

U ( I ) =UL ( I ) 

UOL(I)=U(I) 
soo Continue 

the influence matrix 

IF (INFLU. EQ.O) GO TO 5600 
IF (INFLU, EQ. 2) GO TO 2305 
DO 2300 1=1* NS 

DIZI=DELTc*FLOA'I (IDZ(l) )-DLZ2 
DO 2290 U=l*iJS 

DlZJ=DELTZ*FLOA( (IDZ(J> 1-OLZ2 

K=(J-1)*HS+I 

IREF=I0Z(U) 

ER=ABS(UPG(J)/CUR(J) )*RY 

ER2=Lt<*ER 

IvJ=IABS( 1-U) 

IF (IJ.EO.O.OR.IU.EQ.l XZ> GO TO 2090 
GO TO 2100 

2090 S=ASb(DELTX( IKEFl/ER) 

call sicks*sint*cint> 

EPS1=PI2P*(PIH*(1.0-CCS(S))+C0S(S)*SINT-SIN(S)*CINT) 



Go TO 2250 
2100 continue 

IF (lDZ(I).EU.lUZ(J)) GO TO 2110 
GO TO' 2200 

2110 UlN=«QS(FLOAT(IUX(I)-iOX(J) ) ) 

DlSX=DELTX < iRtiF) *DlN/c.R 
S=Adb(2.0*DISX) 

CALL SICKSrSlNTrClMT) 

EPSl=P I2 P* (SIN ( S ) * (PIh-SInT) -COS (S> *C I NT) 
EpS(K)=DELTX(IR2F)*EPsl/ER 
GO TO 2290 
2200 CONTINUE 

IF <I0X( I) .EO.IQX(J) ) GO TO 2205 
GO TO 2260 
2205 N£P= l 

2210 DIN=a6S(DIZI-DIZJ) 

IF <WEP.L3.2> D1N=DIZI+DIZ0 
disz=beta*din/ek 

OlX=ABS(XW(I)-Xrt(J) ) 

OlSX=DIX/ER 

UX2=oISX*l)ISX 

OISR=50RT(DX2+OISZ*DISZ) 

DS=2-0*ABS(DISR) 

DS2=0S*DS 

CALL CYLIDS.FHOrFHltF^lfFNZ) 

Gl=Pi/DS 

G2=P12-FH1+FiU 

G3=FHo-(FH1+FN1)/DS+FN2 

EFUnC=G1*< (4.0*UX2/OS?-1.0»*G2+<+s 0*QX2*G3/D5) 
EpF=dETA*UELTZ*UELTX ( IREF) *EFUNC/ <Pl4*ER2 ) 

IF (NEP.EQ.2) GO TO 2260 
EPSl=EPF 
2250 NEP=2 

GO TO 2210 
2260 EpS2=EPF 

EpS(K.)=EPSl+EPS2 
GO TO 2290 
2200 EpS ( K } =0 • 0 
2290 Continue 
2300 Continue 

60 TO 2600 

2305 Uo 2500 1 = 1 » iJS 

UlZI=DELTZ*(FLOAT<IDZ{ I) )-DLZ2) 

DO 2^00 U=1»US 

DIZJ=UELTZ*(FLOaT(IDZ( J) ) -DLZ2 ) 

K=( J-l) *NS+I 
IJ=IAaS( I-J) 

Irc.F=I0Z(J) 

DlX=AdS(Xrt(I j-X.,(J) ) 

Ek=Ao5(UPG( J) *RY/CUR) 

ER2=£A*Er< 

DO 2390 JJ=1»4 
EpS 1=0 • 0 

IF (IJ.Eu.O.OK.IU.EQ.lXZ) GO TO 230a 
Go TO 2350 
2300 HEP=1 
2310 AoD=0.0 

IF (Jo.Eu.l) AUj=-DELtX(IREF)/2.0 
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IF ( JJ.Eu.3) ADO=DELTx(IREF)/2.0 

uxsx=idix'+aod)/er 

DX2=0ISX*DISX 

OlZ=AbS(bIZI-DIZJ) 

IF (UEP.EQ.2) OiZ=UIZl+DlZ>J 
AD0=U.0 

IF <JJ .£ 0 . 2 ) AD0=0LZ2 
IF (JJ.Eu.4) aL)U=-DLZ2 
DISZ=3ETA*(DIZ+ADD)/Ek 
0ISR=SQRT(0X2+DISZ*DISZ) • 

DS=2.U*ABS(DISR) 

DSZ=uS*DS 

call cyl(DS,fho.fhi.fni»fn2) 

Gi=Pi/US 

G2-PI 2-FHl+FNl 

G3 =FHu- (FH 1+FN1 ) /05+Fr<2 

EFUNC=G1*( (4.0*DX2/OS2-1»0> *G2+4.0*DX2*G3/DS) 
IF (NEP.EQ.2) GO TO 2360 
EpSl=EFUNC 
2350 NEP=Z 

GO TO 2310 
2360 EpS2=EFUNC 
2390 EP<JJ)=EPS1+EPS* 

EpF=BETA*DELTZ*UELTX(lR£F)/(Pl4*ER2) 

2400 EpS(K)=EPF*(£P(x)+EP(2)+EP(3)+EP(4) )/4.0 

2500 continue 
26oo continue 
Symf=i.o 

IF (ILIFT.EQ.O) SYMF= 2 .0 
Go TO (3000.500Uf 5100) * METHOD 

THE NEwiTON-RAPHbON METHOD 

3000 continue 

3020 Do 3200 1 = 1 * NS 
KKK=l I-1)*NS+I 
SUM=0.0 

DO 3120 J=1 * NS 
K=( J-l ) *NS+ I 
IJ=IABS(I-U) 

IF (1J.EQ.0) GO To 3120 
SUM=SUM+SYMF*EPS(K)*U(J)*U(J> 

3120 CONTINUE 

SUM=SUM+SYMF*EPS(KKK)*ULU)*UL(I) 

AdD=0 .0 

IF (ILIFT.EQ.O) GO TO 3150 
U = I+LXZ 

IF (I.GT.LXZ) I I=I-LXZ 
3140 ADD=0 .5* (U( I ) -UL ( I ) -U( II ) +UL( I I ) > 

3150 FlNTm=UU)-Q.D*um*U(I)-ULlI)+SOM-*-ADD 
KK=N52+I 
A(KK)=FINT(I) 

3200 CONTINUE 

DO 3250 1=1, NS 
DO 3240 J=1 , NS 
K=(J-1)*NS+I 
IJ=IABS(I-J) 

IF ( I J.Eu.O) GO TO 3230 
AU0=0.G 
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IF (IJ.EG.LXZ) ADD=-0.5 
A(K)=2.0*SYMF*U(J)*EP5<K)>ADD 
SO TO 3240 
3230 ADD=0.5 

IF (1LIFT.EQ.0) add=o.o 
A(K)=1.0-U(U)+AOD 
3240 Continue 
3250 continue 

call sauss(a.ns» it) 

Do 3310 I = l»liS 

KK=NS2+I 

U(I)=U(I)-A(KK) 

3310 continue 

NIT=NIT+1 

Do 3320 1 = 1 » MS 

U1 Esl'= ABS COM » -UOL ( I) ) 

IF (UTEST.ST.ACC) GO TO 3330 

3320 continue 

go TO 5800 

3330 continue 

IF (NIT.GE.15) GO TO 5800 
DO 3340 I=1»NS 
U0L(I)=U(I) 

3340 Continue 

Go TO 3020 

The method of parametric differentiation 

5000 PrMT(1)=0.0 
PrM f ( 2 ) -l . 0 
PrMT(3)=0.2 
PrMT( 4)=0.1 
NpD=0 

Do 5010 1=1 # NS 
Y ( I ) =UL ( I ) 

5010 DERY(I)=1.0/FLOAT(NS) 

Call hpcg<prmt.ns,ihlf*aux) 

Do 5020 1=1. NS 
U0L(I)=Y(I) 

5020 U ( I ) =Y ( I ) 

Go TO 5800 

the method of steepest descent 

5100 Do 5050 NZ=1.LZ 
IFG1=0 
IFG2=0 

DC 5040 NX=1»LX 
K= (NZ-1 ) *LX+NX 

kk=k+lxz 

IF (IFG1.EG.1) GO TO 5030 
IF (U(K) .GT.UCR) IFG1=1 
XSU ( NZ ) =XW t K ) +DELTX ( N* > /2 . 0 
5030 continue 

IF ULIFT.EQ.O) Go TO 5040 
IF ( IFG2.EU. 1 ) GO TO 5040 
IF (U(KK) .GT.UCR) IFG2=1 
XSL ( NZ ) -X ti ( KK ) +OELTX ( wZ ) /2 • 0 

5040 continue 



5050 continue 

DO 4d50 NZ = 1 »CZ 

IFG1=0 

IFGZ=0 

DO 4640 NX=lfLX 

fJNAsLX-NX+1 

K;(NZ-1>*LXMMNX 

kk=k+lxz 

IF (IFG1.EQ.1) bO TO 4830 
IF (U(K> .GT.UCR) iFGlrl 
XXU(NZ)=Xw(KJ+0.5*DELTX(NZ) 

4830 CONTINUE 

IF (ILIFT.EQ.O) GO TO 4840 
IF UFG2.E0.1) oO TO 4840 
IF (U(KK) .GT.UCR) IFG2=1 
XXL(NZ)=XW(KK)+0.5*OELTX(NZ) 

4840 CONTINUE 

4850 continue 

WRITE (6* SHOCK) 

C 

Nh=NS/2 

LmIN=1 

lmax=ns 

IMIN=1 

IMAXSNS 

NIT=0 

4000 AC=0.0 

IF (ILIFT.EQ.O) GO TO 4 020 

SUMl— l) « 0 

SUM2-O.0 

DO 4010 1 = 1 * NH 

II=I>NH 

SUM1=SUM1+U ( I ) -UL ( I ) -u ( 1 1 ) HJL ( II ) 

SUM2=SUM2+UL ( I ) *UL ( I ) -UL (ID *UL (ID 
4010 continue 

AC=0.5*SUM1/SUM2 

4020 continue 

Do 4173 I=LMIN»LMAX 

SuM= 0,0 

SuMl=0.0 

KK=(I-D*NS+I 

00 4151 J=lMIf4* I MAX 

K=iu-1)*N5+I 

IF (U.EQ.I) GO TO 4151 

SUM1=SUM1+SYMF*EPS(K)*UI J)*U(J) 

4151 continue 

SuMl=SUMl+SYI'iF*t.PS ( KK ) *UL { I ) *UL ( I ) 

CC=0.0 

IF (ILIFT.EQ.O) GO TO 4152 
1I = DNH 

IF (I.GT.NH) II=I-NH 
CC = U4< I ) *UL ( I ) -OL( I I ) *UL( I I ) 

4152 reK=i,0+2.0*(SU>U-UL<I)+CC*AC) 

TERsSQKT(AUS(TEK) ) 

ZGN=1,0 

IREF=IDZ(I) 

IF (I.GT.UXZ) GO TO 4155 

IF <X<V(D.GT.XSU(IKEF) .MND.Xtfm.LT.XXUtlREFD ZGN=-1.0 
GO TO 4lo0 


34 



4155 CONTINUE 

IF (XW(I).GT.XSL<IREF).AN0.XW(I).LT.XXL(IREF)> ZGNr-1.0 
4160 SuN=2.U*(U(I)-1.0+ZGN*TER) 

SUWM=0.0 

DO 4172 J=LMIN*LMAX 
IF (J.EQ.I) GO TO 4172 
KK=(0-1)*NS+J 
KKK=(I-1)*NS+J 
SuM2=0.0 

DO 41b5 FiJ— 1MIM » IMAX 
Ks (MJ-1 ) *NS+J 
IF (i*lJ.EU. J) GO TO 41 h 5 
SuM2=SUM2 + SYi-tK*EPS(K)*U(MJ)*U(MJ) 

4165 Continue 

SuW2=SUM2+SymF*£PS ( KK ) *UL { J } *UL { o ) 

Cc= 0.0 

IF (ILIFT.Eq.O) GO TO 4166 
JJ=J+NH 

if (u.gt.nh) jj=j-nh 

CC=UL( J)*UL( J)-oL( JJ)»UL( JJ) 

4166 T£K=1.0+2.0*(SUh 2-UL( j)+CC*AC) 

TeR=SQRT(ABS(TEH) ) 

zkn=i.q 

JrEF=IDZ(J) 

If (J.GT.LXZ) go TO 4170 

IF (Xa( J) .GT.XSU< JREF) .AND.XW(J) .LT.XXU(JREF) ) Z«N=-1.0 
Go TO 4171 

4170 Continue 

If iX.V(J) .GT.XSL(JREF) .AND.XW(J) .LT.XXL(JREF) ) ZKN=-1.0 

4171 SU4=2.0*(U< J)-1.0+ZKN*TER) 

Su5=ZKN/(2.U*TEU> 

Su6=4.0*U( I >*SYMF*EPS(KKK) 

SUMM=SUMM+SU4*SU5*SU6 

4172 continue 

4173 OFUdJrSUM+SU^M 
C 

SU^2=0.0 

DO 4169 I=LMII'J»UMAX 
Kk=(1-1)*NS+I 

Su Nl=0.0 

DO 4176 JzlMIN.IMAX 
K=(J-1)*NS+I 
IF (J.EQ.I) GO TO 417 m 
S uMl-5UMl+SYMF*tPS (K)*U(J)*U(J) 

4178 CONTINJE 

SUMl=SUMl + SYMF*EPS<KK> *UL{ I ) *UL( I ) 

CC=0.0 

IF (ILIFT.EQ.O) GO TO 4179 

II=I+NH > 

IF (I.GT.HH) I 1=1— NH 

CC=UL ( I ) *UL ( I ) -UL (II) *UL (II) 

4179 TER=1.0+2.U*(SU i 'H-UL(I)+CC*AC) 

T£K=iQR T ( A6S ( TER ) ) -> 

ZGNsl.O 

IREF=I3Z(I) 

IF (I.bT.LXZ) GO TO 4l&0 

IF (XV( I ) .GT.XSUt IREF) .ANQ.XW(I) .LT.XXU(IREF) ) ZGN=-1.0 
GO TO. 4ldl : ili i • 

4180 continue 
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IF (Xv/(I) .GT.XSL(IREF) .AND.XW(I) .LT.XXL(IREF) ) ZGN=-1.0 
<+181 TEf<l=U< I >-1.0+ZGN*TER 

4189 SuM2=SUF1Z+TER1*TER1 
C 

FFF=SUM2 

SuM=0.G 

DO 4190 I=LMIN.LMAX 

4190 SUM=SUM+DFJ{I)*uFU(I> 

ell=fff/sum 

DO 4200 1SLMIN.LMAX 
U(I)=U(I)-ELL*UFU(I) 

4200 continue . 

nit=nit+i 
D o 4700 1 = 1 *i'iS 
UlEST=A9b(U(I>-U0L(I) ) 

IF ( JTE5T.GT.ACC) GO TO 4710 

4700 continue 
G o TO 5BU0 
4710 continue 

IF (NIT.GE.25) GO TO 5800 
Do 4720 I=1»NS 
4720 U0U1)=U(I) 

Go TO 4000 

WRITE OUT THE solution FOR THE WING 
5800 N*K=1 

IF (NRUN.GT.O) GO TO 5010 
Go Tu (5810*5890*5910)* IWl 
5810 WRITE (6*5820) 

5820 FORMAT (1H0*3X»11H PLaNFORM- ) 

WRITE (6*5830) 

5830 Format (1H .3X.I0H PROFILE- ) 

GO TO 5950 
5890 WRITE (6*5900) 

5900 FoRMaT (1HO,5X*20H PLaNFORM- RECTANGLE ) 

GO TO 5950 
5910 WRITE (6*5920) 

5920 format (iho,3x*i9h planform- triangle > 

5930 WrITE (6*5940) 

5940 FoRMaT ( 1H *3X*23H PROFILE- PARABOLIC ARC ) 

5950 WRITE (6*5960) TAU*SENl*AR 
5960 Format uh *jX* 

1 18H THICKNESS RATIO = F12.5 / 4X* 

2 17H SEMISPAN /CHORD = F13.5 / 4X* 

3 15H ASPEC1 RATIO = Fi5.b) 

6010 nRlTL (6 *6020 AM * BETA * XM * CPST* EMR 

6020 Format <1H0*53X*26H The SOLUTION FOR THE WING // 4X * 

1 15H MACH NUMBER = Fib. 5 / 4Xi 

2 8H BETA = F22.5 / 4X, 

3 22H MACH NUMBER (STAR) = F«.5 /4X* 

4 13H CP (STAR) = F17.5 / 4 X 1 

5 6H XI = F24.5 // 41X, 

6 52H VELOCITY AND PRESSURE DISTRIBUTION ON UPPER SURFACE /) 

6025 WRlTc (6*6030) 

6030 Format uho*4X* 

1 2 H I*5X*8H x/Ch0rD,7a*2H U»11X*4H U+1*7X*7H U ( RED) * 9X * 2HCP* 8X * 

2 8H CP(REO) ,6X*5H P/P(j»10x*2H M*8X*8H M(STAr) /) 

DO 6200 J=1*LZ 
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NN={J-1)*LX+1 

ZS=Zw(NN)/S£mI 

write ( 6 » 6040 ) Zw(NN),ZS 

6040 FoKMrtT ( 1H0 , 3 X » l 9H SPANalSE STATION, Z/CHORD = FS.5,2X, 

1 15H (Z/SEMISPAN = F6.S,1H) /) 

00 6200 1=1, LX 
K=(J-1)*LX*I 

kuk=k+lxz 

IF (NwR.EQ.2) K=KUK 
UwR=U(K) 

Uw=UwR/RU 

Ww=UPG(K) 

WwR=ww*RU 

RTT=Uw*1.0 

WTT=w«+1.0 

Call vecon ( uw » cp , rtm , „ms , ppp ) 

Call vecon(ww»wp»wtm, v< ms#wpp) 

CpR=AMT*CP 

WPR=AMT*wP 

XLCHL=(FLOAT(I)-0.5)/FL0AT(LX) 

WRITE (6,6050) 1,XLCHL#UW,RTT,UWR,CP,CPR,PPP,RTM,RMS 
6050 FORMAT ( 5X » 12 , IX • 9 ( IX , E12 . 5 ) ) 

WRITE (6,6060) i«W,WTT,WWR»WP»WPR»WPP»WTM,wMS 

6060 Format (2ix,a(ix»Ei2.5) ) 

6200 Continue 

NwR=NWR+1 

IF (NWR.GE.3) GO TO 8000 
IF (ILIFT.EQ.O) 60 TO 8000 
WRITE (6,6300) 

6300 Format (// 39X, 

1 53H VELOCITY AND PRESSURE DISTRIBUTION ON LOWER SURFACE /) 

GO TO 6025 
C 

8000 WRITE (6,8030) 

8030 FoRmmT (1H0»48X. 

1 36H THE LINEARIZED SOLUTION IS BASED ON ) 

IF (lwING.GT.O) GO TO 8050 
WRITE (6,8040) 

8040 Format (1H ,50X,32H CARMICHAEL’S NUMERICAL ANALYSIS /) 

GO TO 8070 
8050 WRITE (6,8060) 

8060 Format (ih ,52x,2sh holme • s analytical analysis /) 

8070 continue 

WRITE (6,9000) 

9000 Format (iho»4ox, 

1 4lH THE NONLINEAR SOLUTION WAS OBTAINED WITH) 

Go TO (9010,9030,9050} , METHOD 
9010 WRITE (6,9020) 

9020 FoRMhT ( 1H ,S6X.20H ThE NEwTON»S METHOD) 

GO TO 9070 

9030 WRlTi- (6,9040) uP 

9040 FORMAT (1H ,4bX,4lH Tn£ METHOD of PARAMETRIC DIFFERENTIATION / 
1 55X»19H AT THE VALUE DP = F6.4) 

Go TO 9085 
9050 WRl Tc. (6,9060) 

9060 FORMAT (IH ,47X,31H ThE METHOD OF STEEPEST DESCENT) 

9070 white ( 6 , 90Q0 ) NlT 

9080 FcRM’Al (IH ,56X, 7H AFTER I2,11H ITERATIONS ) 

9085 CONTINUE 



NrUN=NRUN+1 

IF (.4RUN.6T.MOELM) GO TO 9090 
Ak=AFi+D£LM 
Go TO 215 
9090 IwlNG=9 
Go TO 10 
9999 continue 
stop 
end 


subroutine sicks, sint*cinT) 

c 

S2=S*S 

SS=S2 

S4=S2*S2 

S&=S4*S2 

S8=S6*S2 

If (s.gt.i.o) go to 20 

S3=SS*S 

Ss=S3*SS 

S7=S5*SS 

Sg=S7*SS 

S3=S3/18.0 

S5=Sb/600,0 

S7=S7/35280.0 

Sg=S9/32b5920.0 

SINT = S-S3+S5-S7+59 

52 = 52 / 4.0 

S4-S4/96 • 0 

S6=S6/43 20.0 

Se=Sb/322560.0 

ClNJ = O.57721566+ALOG(S>-S2+S4-S&+S8 
GO TO 30 

20 Fi=S6+3a.0272o4*S6+26b.U)7033*S4+335.67732*S2+3S.l02495 
F2=S6+40. 021433*56+332. 62491 1*S4+570.23628*S2+157. 105423 
Fo=Fi/(S*F2) 

6 l=So+42. 242655+56+302. 757 865*S4+352. 01 84 93*S2+21. 821899 
G2=Se+43. 196927*56+462. 465964*54+11 14. 978885+S2+449. 690326 
G0=G1/ (S5*G2) 

SlNT=1.5707963-F0*CO5(S)-G0*5lN(S) 

ClNT=FO*5IN(S)-GO*COS(S} 

30 continue 
retukn 
End 
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subroutine cyL(S»fho#fhi#fni»fn2> 

Pi =3 •1415926535 
PI2=0. 63661977 
IF (S.LT.4.0) GO TO 40 

Ti=4.U/S •>'. 

T2=T1*T1 

PO=( < < (-O.U000037043*T2 + 0.0000173565)*T2-0.0000487613)*T2 
1 +0. 00017343) *T2-0. 001753062 )*T2+0. 3989423 
Q0=( ( ( (0. 000 0032312*T£-0. 0000142073 )*T2+G. 00-00342468) *T2 
1 »0.000086979l)’»T2+0.o004564324)*T2-0. 01246694 
Pl=( l ( ( .U000G42414*T2-. 000 020092 )*T2+. 0000580759) *T2 
1 t. 000225203) *T2 + .002s.21626)*T2+. 3989423 
Gl=( ( ( (-.000003 d 594*T2+. 00 00 1622 )*T2-. 0000398708) »T2 
1 +.000106474l)*T2-.OOu6390400)*T2+. 03740084 
A?2.0/S0RT(S) 

B=A*Tl 

C=S-0. 7853962 
YO=A*PO*5IN(C)+d*QO*CoS(C) 

Yi=-A*Pl*COS(C)+B*Ql*SlN(C) 

GO TO 90 
40 SS=S/2. 

S2=SS*SS 

T=ALOG ( SS )+. 5772156649 

SUM=0 . 0 

TEKM=T 

ro=t 

DO 70 L=1 f 15 
IF (L- 1) 50 » 60 » 50 
50 SUM=SUM+1.0/FLOAT(L-l) 

60 Fl=FuOAT(L) 

T S-T-SUM 

TERM= ( TERM* ( — S2 ) /Ft **2 ) * ( 1 . 0-1 . 0/ ( FL*TS) ) 

70 Y0=Y04TERM 
T£RM=SS*(T-,5) 

SuMri) .0 
Yl=TtRM 
DO 80 L=2»16 

SUMsSUM+l. /FLOAT (L-l) , . 

FL=F40AT(L) 

fli=fl-i. 

TS=T-SUM 

TERM=(TERM*(-S2)/<FU*FL> >*( (TS-.5/FL)/(TS+.5/FLl) ) 

80 Yl=Yl+TERM 

Y0=PI2*YU . 

Yl=-Pl2/5+Pl2*Yl 
90 Y2=2.0*Y1/S-Y0 
FfJlsYl 
Fn2=Y2 

IF (S.GT.11.0) GO TO i20 

SM=S/2.0 

SM2=SH*SH 

AQ = 6 • 0/ ( 3 . 0*Pl ) 

60-4 • 0/P 1 

A(m=-4.0*SH2*AO/15,0 

bN=-4.U*SH2*HO/9.0 

SuMzAO+AN 



ZuM=80+BN 
DO 110 1=1*16 
FLl=f-'LOAT(I) 

AN=-bh2/( (FL1+2.5 )*(FlI+ 1.5) )*AN 
SN=-SH2*bN/(FLI+1.5)**2 
SUM=SUM+AM 
110 ZU^ZU.I+BN 
FhOsSH*ZUM 
Fh1=SH2*SUM 
GO TO 130 
120 S 2 =S*S 
S3=S2*S 
S4=S2*S2 
S5=S4*S 
Sb=S2*S4 
S 7 =S 6 *S 
Sfl=S2*S6 
S 9 =S 6 *S 
SlO=b 2 *S 6 
su=sio*s 
Sl2=b2*SlO 

Fh0=YO+P12*(1.O/S-1.0/S3+9.O/S5- 4 

1 225. 0/S7+1 1025. 0/S9-h93025.0/Sll) 
FHl=FNl+PI2*ll.u+1.0/S2-3.0/S4+45.0/S6-1575.0/S8+ 
1 99225. 0 /SI 0-9823275. 0/S12 ) 

130 continue 
f<£TU KN 
END 


function acsh(Ca) 


cq=ca*ca 

IF (CA.GT.1.0) GO TO 40 

C 3 =Cu*CA 

C5*CO*C3 

C7=C0*C5 

C9=CU*C7 ' 

Cll=CG*C9 

ACSH=CA-C3/6.0+3.0*C5/40.0-15.0*C7/336.0+ 

1 l05.0*C9/345b.u-945.u*Cll/42240.0 
GO TO 50 
40 C4=C0*C0 
C6=Cb*C4 
C 6 =Cu*Cb 
ClO=CO*Cb 

ACSH=ALOG ( 2. *CA) +1.0/( 4. 0*CQ) -3.0/(32.0*04)+ 

1 lb, u/(2bH.0*C6) -105. u/( 3072. 0*C8) +945. 0/( 38400. 0 *C’ 10 ) 
50 continue 
RETURN 
End 
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SjJbROUTInE GAUSS(A»M* IT) 


DIMENSION A ( ,N * 1 ) 

REAL MAX 
IT=0 
£; 1 .£-8 
Nn=N+ 1 
DO 12 I = 1*N 
IF d-1) 3*3*1 

1 DO 2 J=I*N 
M=I-1 

Do 2 K= 1 *M 

2 A(J,I)=A(J*I)-A(J.K)*a<K*I> 

3 KEr=I 

MAX=A3S(A(I»D) 

Do b K=I*N 

IF (MAX-ABS(A(K.I) ) ) 4 »b*5 

4 KEY=K 

MAX=ABS(A(K, I) ) 

5 continue 

IF (MAX-E) 6*7*7 

6 IT=1 
return 

7 Do 8 K= 1 *NN 
maX-a ( i » k ) 

A ( I # K ) =A ( KEY » K ) 

8 A(KEY»K)=.MAX 

ii=i --i 

M=I-1 

DO 12 J=II*NN 
IF (1-1) 9*11*9 

9 DO 10 K=1*M 

10 A(I,J)=A(I*J)-A(IfK)*A(K*J) 

11 A(I*J)=A(1»J)/A(I»I) 

12 continue 

DO 14 I-2*N 

JpN+l-I 

JJ=J+1 

DO 13 K=JU*N 

13 A (JtNN)-A (Jt NN)-A(JtK ) *A (K *NN) 

14 Continue 
return 
End 



I 


subroutine; hpcgiprmt*ndim,ihlf*aux) 

dimension prmkd *aux<16*i> 

common /dy/ o£ry(ioo),y(ioo) *ul(ioo) *dp 

n?i 

IhLF=0 

X=PRMT(1) 

H=PRMT(3) 

PRMT(5)=0. 

DO 1 I=lfNDIM 
AUX(16»I)=0. 

AuX(15#I)=DERY(I) 

1 AuX(l*I)=Y<I) 

IF (H*(PRMT(2)-X) ) 3*2*4 

2 IhLF=12 
GO TO 4 

3 lhLF=13 

4 Call fct(x.ndIM) 

call outp(x»ndim»ihlf) 

IF (PRMT (5) ) 6*5*6 

5 IF (IHLF) 7.7*6 

6 return 

7 DO 8 I=1*NDIM 

8 AUX(8.I)=D£RY(I) 

ISW=1 

GO TO 100 

9 X=X-h 

DO 10 I=1.NDIM 

10 AUX(2.I)=YII) , 

11 IhLF=IHLF+l 
X=X-H 

DO 12 1 = 1 r N0I>-1 

12 AUX(4. I )=AUX(2. 1) 

H=0.5*H 

N;1 
ISW'=2 
GO TO 100 

13 X=X+H 

CALL FCT(X*NDIM) 

N-2 

DO 14 I=1.ND1M 
AUX(2. I )=Y ( I ) 

14 AUX(9. I)=D£RYtI) 

1S*=0 

GO TO 100 

15 DELT=0.0 

DO 16 I=1.NDIM 

16 DELT=DELT + AUX(15* I)*AaS(Y{I)-AUXl4rI) ) 
D£LT = 0 i06666667*DELT 

IF (DELT-PRMT (4) ) 19*19*17 

17 IF (1HLF-10) 11*16*16 
16 1HLFS11 

X=X + H 
GO TO 4 
19 X;X + ri 

call fct(x*nuim) 

00 20 I — 1 » NO I M 
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AUX 1 3* I > =Y C I ) 

20 AuX(10»I)=QERY(i) 

N=3 , 

ISv-=4 

GO TO 100 p 

21 Nsl 

X=X+H 

CaL-L FCT(X.NUlM) 

XrPRMT(l) 

DO 2-i I = i»NDIM 
AuX ( 11 » I ) -OERY ( I ) 

22 Y(l)=AUX(l» I)+H*(0.37b*AUx(&*I)+0.7gib667*AUX(9»I) 

1 -0.20a3333*AUX(10»I)+0.04l66667*DERY(I) ) 

23 XzX+H 
N=N+1 

call fctu.ndim) 

Call outpix,ndik.»ihlF) 

If (PRKT (5) ) 6 » 24 . 6 

24 IF (N-4) 25# 200 » 200 

25 DO 2d I=1»WD1M 
AuX (i'J# I ) =Y ( I ) 

26 AuX(N+7#I)=DERY(I) 

IF (IM-3) 27.29.200 

27 DO 2d I=1.N0IM 
DELT=AUX<9»I)+AUX(9.I) 
delt=delt+delt 

28 Y(I)=AUX(1.I)+0.3333333*H*(AUX(8.I)+DELT+AUX(10»I) ) 

GO TO 23 

29 Do 30 1=1. NO I 
DeLT=AUX(9.I)+AUX(10.I) 

DeLT=DELT<-DELT+DELT 

30 Y ( I ) —AUX ( 1 . 1 ) +0 . 375*H* ( AUX ( 8 » I > +DELT+AUX ( 11 * 1 ) ) 

GO TO 23 

100 Do 101 I= 1 »NQIM 
Z=H*AUX (N+7. I ) 

AuX (5»I ) -Z 

101 Y(I)-AUX(N»I)+0.4*Z 
Z=X+0.4*H 

call fct(z*noim) 

Do 102 I=l»NDlM 
Z=H*JERY(I) 

AUX (6. I )=Z 

102 Y(I)=AUX(N#I)+0. 2969776* AUX(5#I)+0» 158 7596* Z 
Z=X+0.4557372*H 

call fct(z#.nioim) 

Do 103 I — 1 » NU IM 
Z=H*uERY(I) 

AuX(7. i )=Z 

103 Y( I)=AUX(N. I)+0.2ldl004*AUX(5f I)-3.050965*AUX(6. I)+3.832865*Z 
Z=X + H 

call fct(z»ndim) 

DO 104 I — 1 » NjlM 

104 Y( I)=AUX(N» I) +0.l7476u3*AUX(5. I)-0.5514807*AUX(6»I). 

1 +1.205536*AUX(7# I ) +0 . 171 1 S48*DER Y ( I ) 

GO TO (9»13»15»2l)» I GW 

200 IST£P=3 

201 IF.(N-O) 204 . 202 » 204 

202 DO 203 N=2» 7 

DO 203 I=1»NJIM 
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AUX(N-1»I)=AUX(N»I) 

203 AUX (N+o* I )=AUX(N+7* I) 

N = 7 

204 N=N+i 

DO 205 1=1 * NO I M 
AuX(U-l.l)=Y(I) 

205 AuX(N+6*U=DERY(I) 

XzX+H 

IST£P=IST£P+l 
Do 207 I = 1 * ND I M 

D£LT=AUX(N-4» I)+1.333333*H*(AUX(N+6*I)+AUX(N+6»I)-AUX(N+5,I)+ 

1 aUx(N+4*I)+hUX(N+4,I) ) - 

Y(I)=JELT-0.9256193*AUX(16*I) 

207 AuX(16.I)=DELT , 

call fct<x*noim) 

DO 208 1 = 1 * NDlM 

D£LT=0.125*(9.0*AUX(N-1. I > -AUX ( N-3* I > +3. 0*H* (DERY ( I ) +AUX (N+6* I ) ♦ 
1 AUX(N+6*I)-AUX(N+5*I ) ) ) 

AuX(16»I)=AUX(lo* D-DELT 

208 Y(I)=DELT+0.0743817*AUX(16»I> 
delt=o.o 

Do 209 1=1 * NDlM 

209 DELT= JELT+AUX ( 15 * I ) *AsS ( AUX ( 16* I ) ) 

IF (OELT-PRMTU) ) 210,212*212 

210 call fctu*noimj 

call outpu.nuim. ihlfj 

Dp=X 

IF (PRMT(5>) 212*211*?12 

211 IF (IHLF-H) 213*212*212 

212 RETURN . 

213 IF ( H* ( X-PRMT ( 2 ) ) ) 214*212*212 

214 IF (ABS(X~PRMT(2) )-o. 1 *ABS(H) ) 212*201*201 
END 


SuaROUTINE FCT ( X * ND IM ) 

C 

Common /ex/ ex (50) 

Common /dy/ qeryhoo) ,y(iooj *ul(ioo> ,dp 
common /prndtl/ upGduO) * ilift,nrun,symf 
common /hol/ ep( iouooj * a( ioioo) 

IT-9 

Nh=NJIM/2 
Do 1200 I= 1 *NJIm 
S( jM=U . 0 

Do 1120 J=1 * NolM 
K=( J-l ) ♦iJDlM+I 
IJ=IABS( I-J) 

IF (IJ.Eo.O) GO TO UiO 
SUM=3DM+SYMF*EP(K)*Y( J)*Y( J) 

GO TO 1120 

1110 SUM=3U.m + SYMF*EP ( K ) *UL ( J) *UL ( J) 

1120 continue 

k«=nuim*noim+i 

IF (lLIFT.Eo.U) GO TO 1130 
U = I+UH 

IF (I.GT.NH) ii=i-nh 
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A ( KK > =0 . 5* ( Y < I > * Y ( I ) -Y ( I > +UL ( I ) +Y ( 1 1 ) -UL (II)) -SUM 
GO TO 1200 ' ' 

1130 A(KK) = 0»5*Y(I)*Y(I) -SUM 

1200 Continue 

DO 1<250 I = l»f)DI,*l 
Do 1240 U=1»NDIM 
Kr ( J-l ) +NDIM+1 
KK=(U-1)*NQIM+I 
Ij=IABS(I-J) 

IF (IJ.Eu.O) GO TO 1230 
ADOzG.U 

IF ( IA3S( I-J) .EO.NH) a0D=-0.5*X 
AIKK)=2.0*SYMF*a*Y (J> *EP(K)+ADD 
GO TO 1240 

1230 A(KK)=1.0-X*(Y(J)-Q.5) 

1240 continue 
1250 continue 

call gaussu.ndim.it) ■ 

DO 1310 i = l » NU I rt 

kk=noim*nqim+i 

1310 DEHY(I)=A(KK) 
return 
end 


SUBROUTINE OUTP(X.NDIm. IHlF) 

COMMON /DY/ DERY (100) ,Y(100) *UL(100) rDP 

OPTION FOR USER TO PRINT OUT INTERMEDIATE SOLUTION 

RETURN 

EnO 


subroutine vecon(x*cp,am*ams»ph) ■ 

c 

COMMON /VEL/ MACH,MAChSG»BETASQ»BETA»AAM1»'AAM2»AAM3 

real maCH#MACHSG»BETAsQ»BETA 

Xl=X*-i.O 

Xl2=Xl*XI 

CP1-ABS ( 1 • 0+ AAM2* ( 1 . Q-X 12) ) 

Cp=AAMl*(CPl**3.5-l.Q) 

AM=MACH*X 1/SORT (AQS(1.0+AAM2*(1. 0-X 12) ) ) 
AN:S=AV*S0RT(ABS(1. 0/(1, 0 + 0. 16666666*( AM*AM-1. 0) ) ) ) 
PR=AAM3*CP+1,0 
RETURN 

end 



SUBROUTINE carm 

c 

COMMON /VEL/ MACH # MAChSQ # BET ASQ * BET A * A AMI # A AM2 # A AM3 
COMMON /PAR AMS/ N/illNG, PANELS »SREF,REFMOM» CSARt SPAN# OC 
COMMON /WNGPMS/ ROOT<4) »TIP(4) #M#N#TYPE»F(101) »G(10l) # P(101) » 
1 $HEAR ( 101 ) 1 1 SECT #THICK(5) ' "• 

common /HtL/ hhcioo»'1cO) 

REAL MACH # MACHSO » BETA SQ » BETA 

integer panels#oc#type 

logical sym - 

logical thk 

data thk/. false./ 

SyM=. false. 

Data panels/o/ 

Data HALFPl/l.57079633/*PI/3. 14159265/ . 

CaLL if.NGEOM 

n#,ing=panels 

IF (THlCK(l) .GT.0.0.0R.THICK(2) .GT.0.0) ,THK=.TRUE. 

Call eval ithki ' ' 

Nl=NwlNG+l 

Call invert (hh#nwingj 
call force 
panels=o 

return . 

enO 


subroutine wngeom 

c 

Common /vel/ mach,machSQ»betasq#beta»aami> aam2»aam3 
common /params/ nwIng , panels » sref , refmom » cbar » spam » oc 
common /wngpms/ rootu) *tiP(4) »m»n»type#F(1oi> #G(i'oi) »p(ioi) » 

1 sHEARUOI) #ISECT,THIcK(2) »perchd#zt»dzdx# curx . 1 

Common /psings/ pw<ioo) *alphat(ioo> 

common /poata/ xbaR(ioO) »aRea(ioo) *costhuoo) »sinth(ioo).*sym(ioo) 
common /scrap/ x(ioo»4> »yiioo»2) *z<ioo#2) #xcpt<ioi) # 

1 SLOPE ( lul ) #CHRoOT( 202) • ChTIP ( 202 ) » ZUROOT ( 101 > »ZUTlP<i01,> » 

2 ZLROOKIOI) #ZLTIP(101) rSCF(lOl) 

common /curv/ cuR(ioo) , . 

real mach #machsq»betasq» beta 
integer panels# panmax, oc# type 
logical sym, fin 

data panmax /iou/ ‘ . . 

heLl=o. 000001 

IF l M • LE . 0 . OR . N . LE « 0 ) RETURN 

MlSM+1 

Nl=N+l 

F(1)=K00T(1) 

G(1)=TIP<1> 

F(Nl)=ROOT (2) 

G(N1)=HP(2) 

P(1)=R00T(3) i 

P(M1)=TIP(3) 

SHEAR (1)=R00T(4) 
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ShCAK(M1)=TIP(4) 

0|<-RuOT ( 2 ) -ROOT ( 1 ) 

DT=TiP(2)-TIPll) 

SsPAN=T IP ( 3 ) -ROOT { 3 ) 

PAHEA= ( DR+DT ) *SSPAN 

sci=i. 

S C2=l. 

IF (TYPE.EQ.i) GO TO 20 
IF ( I YPE.EG.2) 60 TO 10 
Call fill (f,n> 
call fill t g . n > 

10 IF (TYPE.EU.3) GO TO 20 
call fill (p,m> 
call fill (shear.m) 

20 Xx-TIP ( 3) -ROOT (3) 

00 30 I=l»Nl 
ChR00T(I)=0. 

ChTIP(I)=0. 

SlOPE(I)=(G(I)-F(I) )/ xX 
30 XCPT(I)=F(I)-ROOT(3 )*sLOPE(I) 

DO 50 1 = 1 .N 

Ir < Abf. ( OR ) . L T . hELL ) (0 To 40 
ChROOT( l) = ( (F< 1) <-F(I + !) )/2.-F(l) )/DR 
40 IF (m8S(DT) .LT.hElL) GO TO 50 

ChTIP(I)=( (G( i)+G(I+D )/2.-G(l) )/DT 
50 continue 

IF (UT.GT.O.) GO TO 70 
00 60 1=1. N 
60 ChTIP(I)=CNR0OT(I) 

70 Do 90 J=1 . M 
Y l=P ( J) 

Y2=P(J+1) 

YPER=( (Yl+Y2)/2.-P(l) )/(TlPl3)-R00T(3) ) 

Y ptRl = l • O-YPER 

ThCK=THICM2)*YPER+THICK(1>*YPER1 
DO 90 1 = 1 » N 

panels=panels+i 

If (PANELS. GT. PANMAX) GO TO 210 
X(PANELS.1)=XCPT< I)+Yi*SLOPE(I) 
X(PANELS. 2 )=xCPr(IJ+Y 2 *SL 0 PE(I) 

X (PANELS. 3) =XCPT ( 1+1 ) +Y1*SL0PE ( 1+1) 

x (Panels. 4 ) =xcpt < i+d+y2*slope( i+i) 

Y(Panels.d=yi 

Y(PANELS»2)=Y2 
Z( PANELS. 1)=SHEAR(U) 

Z ( PANELS . 2 ) =S' 'EAR ( J+l ) 

IF (AdS(THCK) .LT.HELL) GO TO 80 
P£RCml)=CNT Ip ( I ) *YPER + CHKOOT ( I ) *YPERl 

call sectin 
C|jR(PAIiElS)=THCi\*CURX 
80 ALPHAT (PANELS) =THCK*DZDX 
FlN=. FALSE. 

If (ADS(Yi) .LT.HELL. Ar,D.ABS(Y2) .LT.HELL) FIN=.TRUE. 
90 SYN(PANELS)=.N0T.FIN 

calculations of wing profile 

KTMAX = Tr(lCK(l)*jR 
riPMAX = T(iICK(2)*DT 
DO 130 1=1. N1 

IF (AdS(KTiiAX). LT.HELL. OR. ABS(DR). LT.HELL) GO TO 100 



oooooooooo 


P£RCH0=CF(I)-F(1) )/or 
' call sect in 

2UR00K I)=R00T(h)+2T*kTMAx 
2LR00K I)=R0uT(4)-ZT*rTMAX 
Go TO 110 

100 2U*00T{I)=K00T(4) 

2lROOT(I)=KOOT(4> 

110 IF (AfaS(TIPMAX) .LT.HE lL.OR.ARS(DT).LT.HELL) GO TO 120 
PERChO=(G( I )-Gd) j/OT 
Call sectin 

2UVIP(I)=TIP(4)+2T*TIf>MAX 
2LTIP(I)=TIP(4)-2T*TIPMAX 
Go TO 130 

120 ZUTlPl I )=TIP(4) 

2LTIP(I)=T1P(4) 

130 Continue 
210 return 
End 


subroutine fill (afill.nfill) 

c 

real AFILL(l) 

IF (NFILL.LE.l) RETURN 

UELi ( AFILL ( NFILL+1 ) -AFILL <1)1 /FLOAT (NFILL ) 
DO 1U I=2.NF1LL 
10 AFILL(I)=AFILL<I-D+DEL 
return 
end 


subroutine sectin 

c 

common /wngpms/ root< 4) * tip<4) #m#n» type.f(iqi) »G(iod ,p(iod . 

1 SHEmRUOI) . ISECT.THICM2) . X . Z . DZDX , CURX 
DIMENSION XTAB ( 19) » ZTABH19). ZPTA8K19). CUTaB1(19) 

the four data Tables given below are for a rae 101 profile 
and they represent 19 chordwIse values for the profile 

X/CHORD STATION TOGETHER WITH THE THICKNESS » SLOPE AND 
curvature measure respectively as given for a profile 
thickness ratio of out. for arbitrary profile input the 
appropriate data table- s must be supplied by the user. 

IF THE NUMbEK OF CHORuKiSE STATIONS USED IS DIFFERENT 
FROM 19. THE DIMENSION STATEMENT ABOVE AND THE THREE 
SEOUtNTIAL CALLiNG STATEMENTS STARTING WITH STATEMENT 100 . 
must also be changed accordingly. 

Data XT Au /O. ». 005. .U125*. 0250. .0S».075». If .15». 2^.25*. 3».4f 
1 • 5 » . fa . » 7 * . d » .9 . .95 . 1 . / 

DAlAZTABi /O. . .0871 . . 1 369. • 19 17 » .2659 . .3191. ,3607» .4220* .4630. 
1. 4dbb.. 4997. .4001, .42u7» *3531 . . 2681 1 789. . 0894 . . 0447. 0 .'/ 

data zptabi /39. 05, e.oas. 5. 205.3. 709.2. 452.1. 863. 1.434.1,001. 

1.6575. .3b75. .0640.-, 393e.-.b50ti'-. B058»“. 8823 .-.6943. -.8943. 
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2-. 6945#-. 8944/ 

Data CuTABI /-7000.0#-420.0»-120.0»-83.1»-32.6»-19.60l »-i2.35# 

1 —8. 102 * -6. 915# -B. 92 > -4 . 9fib# -3.384# r2 .117# “1 .173# -0. 602 » 

2 -0.331»-0.182#-0. 150, -0.125 / 

IF (X.LT.O. .OR.X.GT.l. ) GO TO 130' 

IF (ISlCT.E&.O) GO TO 100 

GO TO ( 10 # 30 # 50 } r 1S£CT , • " ' 

10 z= 2 .*x*(i.-x) 'V 

DZDX=2.*(1.-2.*X) 

CuHX=-4.0 

RETURN 

30 IF (X.GT..5) GO TO 40 . , ‘ .? 

Z=X 

0ZDX=1. 

Cu«X=-4.0 ! “ 

RETURN 1 

40 Z=l.-X 
OzDX=-l. 

CuRx=-4.0 

RETURN 

so continue 

IF (X.GT.0.0) GO TO 6y 

Z=0.0 

DZOX=100.0 , . 

CuRX=-5.0 

RETURN 

60 Sx=SURT(X) 

X2=X*X , 

Z=5.0*(0.2969*SX-0.12e>*X-0.3516*X2+ 

1 0.2b43*X2*X-0.10l5*X2*X2) 

DzOX=5.0^(0.14845/SX-u. 126-0. 7032*X+ ' ! . 

1 0 • 8529*X2-U «40b*X2*X) ' ; ‘ 

CuRX=-5.0*(0.074225/(5X*X)+0.7032- ’ 

1 1.705Q*X+1.21B*X2) 

RETURN 

ioo call taint (xta8»ztabi#x»z#i9»2> 

Call taint (xtau»zptabi»x,dzdx# 19.2) 
call taint (xtau»cutadi»x#curx#19»2) 

return ' 

130 2 = 0 . 

DzOx=0. 

cuRx=o.o 

return ■ . • ■ . 



subroutine taint (XTAb*FTAB*x*FX,N*K) 

REAL XTAb(l) * FTAB( 1 ) , C ( 10 ) * T ( 10) 

DO 10 1—1 , N 

IF (X.GT.XTAB(I) > GO TO 10 
Jzl 

&0 TO 20 

io continue 

J=N 

20 J=J-(X+l)/2 

IF (J.LE.G) J=1 
30 Mr J+K 

IF IM.LE.N) GO TO 40 
J=J-1 
GO TO 30 
40 KP1=K+1 
GO TO 50 

entry tnt 
50 Do 60 L= 1 .KP 1 
C(L)=X-XTA6< J) 

T<L)=FTAB<J) 

60 JrJ+1 

DO 80 J=1*K 
Ir J+ 1 

70 T(1') = (C(J)*T(X)-C(1)*T(J) )/(C(J)-C(X>) 
1 = 1+1 

IF (I.LE.KPl) GO TO 7o 
80 CONTINUE 
Fx=T(KPl> 

RETURN 

end 


SUBROUTINE EVAL (THICK) 

c 

COMMON /VEL/ f4ACH,MACHSU,3ETASQ,eETA*AAM;, AAM2.AAM3 
COMMON /PAR AMS/ Nw ING , PANELS , SREF , REFMOM , C3AR * SPAN* OC 

common /poata/ xbars( ioo) . areauoo) *cosths(iooj *sinths<ioo) , 

1 STMllOO) 

common /comps/ xpr i me , ypr i me * zprime , u * v . v< * b * bterm * eps * sub * bpos * xpm 
it 

common /scrap/ x(ioo*4) *y ( 100*2) *zdoo*2) *xcs(ioo) * ycs(ioo) * 

1 zCS(lOO) *A(100) *XuAR{ 100 ) *XC(100) *YC(100) *ZC(100) * 

2 SlNTH(lOO) *COSTH(100) *UWT(100) 

Common /hel/ HH(ioo*i 1( o) 

COMMON /HOL/ HHi(100*i00) *HH2(100*100) 

CoMt'iON/XZ*/XV(100) *ZV(1U0) 

REAL MACH*i;ACrtSu*BETAsQ*BETA 
REAL LE 

Integer panels.uc 

Logical sub, bpos *bneg, bineg*B2neg*diag* wing* thick *twing 
logical sym 

Data PIl/.3l£33099/*PIp/.159l549/ 

DATA Pl4/7.957747163E-2/*Pl6/3.978b73581E-2Z 
UsELrrO.U 
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ST0TAL=0. 

8NEG=. FALSE. 

DO 10 1=1 *PAuELS 
CR=X(I*3)-X(I»1) 
DELY=Y(I*2)-Y(I»1) 
□EL2=Z(I*2)-2<I*1) 

SPN=SORT l DELY*0 £LY+DEl2*DEL2 ) 
CT-a( I *4)-X( 1*2) 

ARt'A< I) = (CH+CT) *SPN/2. 
STOTAL=ST01 AL+AUEA ( I ) 
YPtR=ll.+CT/(CR+CT) )/3. 
YPER1=1.-YPER 

LE=YPER*X (I#2)+YPERl*x(I»l) 
TE=YPE«*X(I»4)4-YPER1 #x(I» 3> 
XCS(I)=.y5*TE+.05*LE 
XC(I)=XCS(I)/DErA 

XBARS(I)=(LE+TE)/2. 

XV(I)=XBARS(I) 

XBAR( I )=XBARS( I )/BETA 

YCS(I)=Y (I*1)*YPER1+Y(I*2)*YPER 

Zv(I)=YCS(I) 

YC<I)=YCS(I) 

2CS ( I ) =Z ( 1 * 1 ) *YPER1+Z (1*2 ) *YPER 
2C< I)=ZCS(I) 

SINTHS( I)=OELZ/SPN 
SINTH(1)=SIUTHS(I) 

C0STH5( I)=D£LY/^PN 
10 C0BTH(I)=C05THS(I) 
STCIAL=STGTAL*2. 

IF (SREF.LT.O.) SREF=sTOTAL 
IF (CBAR.LT.O.) C8AR=sURT( 5REF) 
EpS=AMAXl (SPAWf CBAR)/iOOOO. 

IF (BETASO) 40 r 40 t 30 
30 SUB= . TRUE « 

X0TEKM=1. 

UC0N=PI8 
UTC0N=PI2/BETA 
VnC0N=BETA*PI3 
VwTCUN=PI 2 
GO TO 50 
40 SUb=. FALSE. 

XBTEKM=-i. 

■ UC0N=PI4 

UTCOil=PIl/OETA 

VwC0(N.=6ETA*PI4 

VwTCO,N1=PI1 

50 CO 2fc>G I=1 »PANEuS 
CoST=COSTH(I ) 

SlNT=bINTH( 1) 

wing=.true. 

TlvlNG=THICK. AiiD.WING 
Xl = X (1*1) /BET A 
X2=X ( I * 2 ) /BETA 
X3=X ( I * 3 ) /BETA 
X4=X< I i 4 ) /BETA 

ri=r(i*n 

Y 2 = Y ( 1*2) 

Zl=2(I*l) 

22=2(1*2) 



0ELTAY=(Y2-Y1)*C0ST+<Z2-Z1>*SINT 

8l=(X2-Xl)/0£LTAY 

B2=(x4-XJ)/D£LThY 

l3lN£o=31.LT.O. 

B2N£b=B2.LT.O. 

Bi=Ad5(31) 

d2=AUS(82) 

BTERi'il = SGRTUaS(3l*Bl + XdTERM) ) 

6TtRf'i2=SuRT ( AdS ( 02*02+ Xb TERM ) ) 

DO 2dU J=l* PANELS 
DIAG=I.E&. J. AND. WING 

Xw=SINT*COSTH( J> , 

Xx=COST*SINTH<J) 

XY=COST*CGSTHl J) 

XZ=SINT*SINTH( J) 

SINTK=XW-XX 

CoSTR=XY+XZ 

SlNTL=Xx+XX 

CoSTL=XY-XZ 

C BEGIN THE EIGHT SEPARaTE SET-UPS FOR COMP AND TCOMP 

BPOS=.NOT.B1NlG 
B=B1 

bterm=btlrmi 
C SETUP 1 

DELTAY=YC(J)-Y1 

DELTmZ=ZC(J)-ZI 

XPRImE=XC(J)-X1 

YpRJ'>iE=DELTAY*C0ST+DELTAZ*SINT 

YP3::YPRIME 

IF (B1NEG) YPRIME=-YPrIME 

zprine=olltaz*cost-oeltay*sint 

ZPJ=ZPHIME 

call comp 

aavr=v 

aawr=w 

IF (DIAG) GO TO 60 
UU=U 

GO TO 70 
60 UU=0. 

70 IF (.NOT. TWINS) GO TO 90 
XPMT=XbAR(U)-Xi 

call tcomp 

UT=U 

IF (DIAG) GO TO 80 
AAVRT=V 
AAwRl 
GO TO 90 
80 AAVRT=Q. 

aawri=o. 

A(jN=V 

c setup 2 

90 DELTAY=-YC( J)-YI 

YPRIME=DELTAY*COST+DElTAZ*SINT 

YP3N=YPRIME 

IF (BINES) YPRIMEs-YPrIME 

ZPRIME=DlLTAZ*CoST-DElTAY*SINT 

ZP3N=ZPRIME 

call comp 
aavl=v 
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Aa*L=W 

uu=uu+u 

IF ( .NOT. TWINS) GO TO 100 

call tcomp 

UT=ur+u 

aavlt=v 

AAftLT^W 
C SETUP 3 

100 DELTAY=YC(J)-Y2 
DELTA2=ZC( J)-Z2 
XPK1,*!E = XC ( J)-X2 

YPRIME=OELTAY*COST+DElTAZ*SINT 

YP4=YPRIME 

IF (BINES) YPRI.-.E=-YPkIME 

ZPRI-ME=DELTAZ*COST-DElTAY*SINT 

ZP4=ZPRIN'.E 

Call comp 

aavr=aavr-v 

AanR=AAWR-W 

IF (l)I AG) GO TO 110 

UU=UU-U 

110 IF (.NOT.TWING) GO TO 130 
XPMT=X8AR( J)-X2 

Call tcomp 
ut=ut-u 

IF (01 AG) GO TO 120 
AAVRT = AAy/RT-V 
AA^RT=AAwRT-W 
Go TO 130 
120 A6N=A3N-V 
C Setup 4 

130 DELTAY=-YC(J)-Y2 

YpRIIME=DElTAY*COST+OElTAZ*SINT 

YP4N=YPRI>iE 

IF (dlNEG) YPRIMEs-YPrIME 

ZpRIME=OELTAZ*COST-OELTAY*SINT 

2P4N=ZPRIME 

call comp 

aavl=aavl-v 

aawl=aawl-w 

UU=Uu-U 

IF < .NOT. TWINS) GO TO 140 

call tcomp 

ut=ut-u 

aavli=aavlt-v 

AawLT=AA,vLT-W 
140 Bp0S=.N0T.B2Ut£G 
B=B2 

BlERM=BTERM2 

c setup 5 

XPKli‘i£=XL( J)-X3 
YpRlME=YP3 

IF (oZNEG) YPRI,'5 E=-YPmIME 

ZPRIME=ZP3 

Call comp 

aavr=aavr-v 

aawr=aa«k-w 

IF (01 AG) GO TO 150 

UU=UU-U 
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150 IF (.NOT.TWING) Go TO 170 
XFMT=XbAR{ J)-X3 
Call tcomp 
ut=ut-u 

IF (OIAG) GO TO 160 
AaVRT=AAVRT-V 
AAWRT = AAwRT-W 
Go TO 170 
160 ABN=ASN-V 

c setup & 

170 YpRlME=YP3N 

IF (U2NEG) TPRIi-iEs-YPrIME 
ZPRIM£=ZP3N 

call comp 
aavl=aa\/l - v 
aawl=aawl-w 
uu=uu-u 

IF (.NOT.TWING) go TO 180 

call tcomp 

UT=UT-U 
AAVLT=AAVLT-V 
AaWLT=AA*LT-W 
C SETUP 7 
180 YPRIME=YP4 

XpRlMe=XC( J)-X4 

IF (U2NEG) YPKIME=-YPkIME 

ZPRIMc=ZP4 

CA'-L COMP 

AA VR=AAVR+V 

AAWR=AAWK+W 

IF < Ul AG) GO TO 190 

UU=UU+U 

190 IF (.NOT.TWING) GO TO 210 
XPMT=XaAR(J)-X4 

call tcomp 

UT=UT+U 

IF (Ul AG) GO TO 200 
AAVRTzAAVRT+V 
AaWRT=AAWRT+W 
Go TO 210 
200 A6N=*BN+V 

c setup a 

210 YPRINE=YP4N 

IF (bzNEG) YPRlwE=-YPnIME 
ZPKIME=ZP4N 

call comp 
aavl=aavl+v 
aawl=aawl+w 
uu=uu+u 

IF (.NOT.TWING) GO TO 220 

call tcomp 
ut=ut+u 
AAVLT=AAVLT+V 
aawlt =aa*lt +w 
220 CONTINUE 

A( J)=(AAVR*SINTR+AAVL*SINTL+AAwR*COSTR+AAWL*COSTL)*VWCON 
Uv.PVAL = UU»UCON 
IF (.NOT.TwING) GO TO 240 

uwT(U)=ut*utcon 



Go TO 245 
240 U*T(J)=0. 

245 Hh(J»I)=A(J) 

HhI ( i*J) = UwPVAL 
Hh2(I»J)=UWT(J) 

250 continue 
2&0 continue 

IF (.NOT. THICK) WRITE <6>370) 
return 

370 FoRMhT (49H0WING THICKNESS MATRIX CALCULATION WAS SUPPRESSED) 
end 


subroutine comp 

c 

common /comps/ a*y#z.u*v#w»b»bterm»eps.sub»bpos>xt 
logical suu.aPos 

data pi /3. 14159265 / .haLFPi/i. 57079633 / 

Data p 132 / 4 . 712389 / »zlro/o. a owe/ 1 ./ 

D2=0.0 

10 Continue 

IF (A8S(B) .LT.1.0E-S) GO TO 150 

X2=X*X 

r 2 =Y*r 

IF (ABS(Z) .LT.EPS) GO TO 120 

22= Z*Z 

R2=Y2+Z2 

Br2=B*R2 

D=SukT(X2+R2) 

A=X-t)*Y 

RpRIM£=SuRT(A*A+BTERM*BTERM*Z2) 

Fb=(X+U)/R2 

F 2 =ALOG(AbS( (B*X+Y+BT£RM*D)/RPRIME) ) 

U=AT*N2 ( 2*U » QR2-X* Y) + aT AH ( Y/Z) 

V;Z*F 6— B*U 

W=BT£RM+F2-Y*F6-B*ALOo(ABS< (X+D) *RPRIME/BR2) ) 

Go TO 190 

120 A=X-o*Y ' 

AA=AoS(A) 

D=SQRT(X2+Y2) 

F2=0.0 

. IF (AbS(d*X4-y+orERM*0),LT.1.0E-8) GO TO 121 
F2=AL0GU8S( (b*x+Y+UT£RM*0)/AA) ) 

121 Continue 
u=halfpi 

IF ( A . LT . ZERO ) GO TO 140 
IF (Y.GT.ZtRO) GO TO 130 
Uz-HALFPI 
GO TU 140 
130 U=PI32 
140 V=-B*U 

W=DTbRyi*F2-(X+D)/Y-B*»L0G(ABS( ( X+ D > * AA/B/Y2 ) ) 

GO TO 190 

150 IF (ABS(Z) .LT.EPS) GO TO 160 
X2=X ♦ X 
Y2=Y*Y 
Z2=Z*Z 
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R2=Y2+Z2 
DrSQHT (X2+R2) 

U=A1AN2(Z*u»-X*Y >+ATA,-j(Y/Z) 

F6=(X+u)/R2 

V=Z*Fb 

W=ALOG(AdS( (Y+D)/SQRT(X2+Z2) ) )-Y*F 6 
GO Tu 190 

160 OzSGiKT (X*X+Y*Y) 
u=halfpi 

IF (X.LT.ZERO) GO TO 18 O ' 

IF ( Y .G r .ZERO) GO TO i7Q 
U=-HALFPI 
Go TO 190 
170 U=Pli2 
180 V=ZEKO 

W=ALOG(AdS( (Y+D)/ABS(X) ) )-<X+D)/Y 
185 continue 

190 IF (dPOS) RETURN 
U;-U 
w--w 
RETURN 

entry tcomp 

X 2 -XT *XT 

Y 2 =Y*Y 

A=XT-B*Y 

IF (ABS(Z) .LT.EPS) GO TO 270 

Z2-Z*Z 

R2=YZ+Z2 

RrSOR T ( R2 ) 

OrSuRT ( X2+R2 ) 

RpRIME=SGRT(A*A+BTERM*6TERM*Z2) 

F2=AL0G(ABS{ (d*AT+Y+BTERM*D)/RPRlME) )/BTERM 

U--F2 

V/ = d*F2-AL0G(ABS( (XT+O) /R ) ) 
v*=ATAN2(Z*D,B*R2-XT*Y) , 

RETURN 

270 D=SQRT(X2+Y2) 

F2=0«0 

IF (AdS(d*XT+Y+tiTERM*u) .LT.1.0E-8) GO TO ,271 
F 2 =ALOG(ABS( (B*xT+Y+dTERM*D)/ABS(A) ) ) /BTER'^ 

271 continue » . . 

W=ZERO 

IF (A*Y.GT.ZERO) w=pi ' 

U=-F2 

V=B*F2-AL0G(A6S( (XT+D) /ABS(Y) ) ) 

RETURN 

END 
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SUBROUTINE INVERT (AA.NN) 

real aauoo.ioo) 
nnn=nn 

Do 30 1 = 1 , NN 

ITEST=I 
PI VOT=AA ( I , I ) 

AA ( I » I ) =1 • 

DO 10 L=1 * NN 
LTfc-Sl =L 
ATEST=AA(I#L) 

10 Aa(I ,L)=AA( l.LJ/PIVOT 
DO 30 M=1»NN 
MjES T =M 

IF (W.EQ.I) GO TO 30 
TT=AA(M»I) 

AACM* 1)=0. 

DO 20 L=1 ,NN 

20 AA(M»L)=AA(M»U-AA(I»L)*TT 

30 continue 
return 
end 


subroutine force 

this subroutine includes the namelist /lift/ 

WHICH SPECIFIES THE GlGBAL ANGLE OF ATTACK (IN DEGREES) 

OF THE WING PLANFOKM aND THE CAMBER OR LOCAL ANGLE 
OF ATTACK (IN RADIANS) OF EACH WlNG PANEL. 

common /v el/ mach , machSO , bet asq * bet a » aami , aam2 * a am3 
Common /H ARAMS/ Nw I NG, PANELS .SREF.REFMOM.CBAR, SPAN » 0C 
COMMON /PSlNGS/ Prf( lot) »ALPHAT(100> 

common /pdata/ xbAR(iuO) .areauoo) .costhuoo) »sinth( ioo) »symi ioo) 

COMMON /SCRAP/ UW(IOO) , DUMMY (2001 »CAMBER(100) * ALPHA ( 1 00 ) * CP ( 100 ) , 
1 AdOO) rLt(lOO) ,CPU(10u) rCpL(lOO) »DELCP(1UO)*UWT(100) 

common /hel/ hhuooYuo) 

common /hol/ hhi(ioo» ioo) »hh2(ioo,igo) 

common /prnqtl/ upg(i 0 o> , ilift.nrun.symf 

REAL MACHfMACHSO*BETASQ*BETA 

Integer panels. oc 

Data R A D/57, 2957795/.* LF A/0.0/ 

namelist /lift/ ilift, camber, alfa. 

IF (NRUN.GT.O) GO TO 50 
read (5, lift) 

WRlTt. (6.LIFT) 

50 AnGLE= ALFA /RAD 
ARADEGzALFA 
WRITE (6,200) ARAOEG 
IF (NvINO.Ea.O) GO TO 60 
DO 60 1=1 » NWING 
60 Pw(I) =0 « 

DO 70 J=l, NWING 
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AxX=-C amber ( J ) -angle*costh ( J ) 

Do 70 I=1»NWING 
70 Pk,(I)=Pft'(I)+HH(I.J)*AxX 
ao Continue 

do 90 i=i*panels 

UwT(I)=0.0 
90 U W (I)=0.0 

Do 1 UU J =1 »N,vINo 

PwTj=ALPIiAT(J) - 

ClJ=Pw(J) 

Do luO .1=1 *PANElS 
UwT(X)=UWT(I)+HH 2(JH)*PWTU i 
Uto ( I ) =UW ( I ) +HH1 { J, I ) *CLJ 

ioo continue 

Do 140 1 = 1 . NW ING 
Ux=UW(I)+UWT(I)+Pw(I)/4. 

UpG ( i ) =UX 

IF (ILIFT.EQ.O) GO TO 140 
Ux=Urf(I)+UwT(I)-PW(I)/4. 

I I = I+N«) InG 
UpG (II) =UX 
140 Continue 

RETURN 

200 Format (4X»3lH CONFIGURATION ANGLE of ATTACK=»F7.3»8H DEGREES) 
End 



APPENDIX C - SAMPLE INPUT AND PRINTOUT 


1 ■ - •• If. 

The subject computer program is dimensioned for a total of 100 panels to represent the 
wing planform and this limits the number of panels which can be used for the lifting case to 
50. A symmetrical nonlifting wing, however, can be represented by a maximum of TOO 
panels. The storage requirement is approximately 145,000 in octal . The input data which 
defines the wing geometry and its relative position to the oncoming flow of given strength 
are arranged in the following sequence of namelists. These namelists are read in from the 
main program except the namelist LIFT which is read in from subroutine FORCE. . 

Note: The given numbers in the input data correspond to the sample calculation case 
shown in the printout. Namelist $ SHOCK is only used for the output of the 
calculated locations of the shock and sonic lines in connection with METHOD = 3. 

$ WING 

IWING = 0 Type of wing, i.e., arbitrary input (0), rectangular planform with 

parabolic arc cross-section (1), triangular planform with parabolic 
arc cross-section (2). 

Note: An arbitrary half-wing planform is defined as a trapezoidal 
planform with the two parallel sides in the direction of the axis 
of symmetry (i.e., the x-axis). The coordinates of the four 
corner points is the required input, see ROOT(i) and TIP(i), 
together with the two thickness ratios THICK(l) and THICK(2). 

INFLU = 1 Type of influence matrix to be used, i.e., semi two-directional 

(1), complete planar (2). INFLU = 0 for linearized solution only. 

Note: It is recommended to use the option INFLU = 1 unless the 
program has been modified to account for a large number of wing 
panels, i.e., N » 200-300. 


METHOD = 1 Method of calculation, i.e., Newton (1), parametric differentia- 

tion (2), method of steepest descent (3). 

Note: METHOD = 1 should be preferred for purely subsonic flows 
whereas METHOD = 2 applies to slightly supercritical flows. For 
the calculation of discontinuous flow the option METHOD = 3 is 
the only one applicable. 

ISECT ~ 3 Type of wing section, i.e., arbitrary input (0), parabolic arc (1), 

double wedge (2), NACA 00XX (3). 

Note: The guidelines for arbitrary section input is described in 
Appendix B under subprogram name SECTIN . 

LX = 10 Number of panel division in the chordwise direction. 


59 



LZ = 5 

ROOT (1) =0.0 
ROOT (2) = 1 .0 
ROOT (3) =0.0 
ROOT (4) = 0.0 
TIP (1) = 0.0 
TIP (2) = 1 .0 
TIP (3) = 3.5 
TIP (4) = 0.0 
THICK (1) =0.12 
THICK (2) =0.12 
SPAN = N.A. 
TAU = N.A. 

MACH = 0.72 
DELM = 0.0 
NDELM = 0 


IUFT = 0 

CAMBER = 100 * 0. 


Number of panel division in the spanwise direction 
x-coordinate of the leading edge at the root, 
x-coordinate of the trailing edge at the root, 
z-coordinate of the root section, 
y-coordinate of the root section, 
x-coordinate of the leading edge at the tip. 
x-coordinate of the trailing edge at the tip. 
z-coordinate of the tip section, 
y-coordinate of the tip section. 

Thickness ratio of the root section . 

Thickness ratio of the tip section . 

Span to root-chord ratio for the cases of IWING = 1,2. 
Thickness ratio for the cases of IWING = 1,2. 


Not 

required 

for 

IWING = 1,2 


$ FLOW 

Initial freestream Mach number (M ). 

00 

Desired increment of the freestream Mach number. 

Number of increments desired of DELM. 

Note: The options DELM and NDELM are introduced for the 
convenience of the user to reduce the number of input data cards 
for the case where only a change in Mach number is desired . 


$ LIFT (required only for IWING = 0) 

Type of flow considered i.e., none ire ulatory (0), circulatory (1) . 

Local slope of each planar panel element which defines a cambered 
wing (in radians). 

Note: The subscripted quantity CAMBER (i) is read in as an array 
of index i = (m - 1)LX + m where 1 s.n £ LZ and 1 £ m £ LX and 
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the indexing starts with the panel located at the leading edge of 
the root section . The input can be simplified by use of the notation 
shown for the sample cose. 

ALFA = 0.0 Angle of attack (in degrees). 


The printout lists these namelists for each calculation case. The results from a flow 
calculation are written Out in two rows for each chordwise calculation point. The upper 
row refers to the final nonlinear solution, whereas the lower row refers to the linearized or 
Prandtl-Glauert solution. The calculation points are located at the centeroid of each wing 
panel and results are given for each spanwise station in the following notation: 

| Chordwise sequential number of wing panel 


X/CHORD 

U 

U + 1 
U (RED) 

CP 

CP (RED) 
P/PO 

M 

M (STAR) 


Chordwise location of panel centeroid with reference to local 
chord length 


u - u 


Normalized perturbation velocity, = 


Normalized velocity, - 


Reduced perturbation velocity. 


(x + 1)M u - u 


Pressure coefficient, ^ 

xM 


1 - M 

00 

1 + l (h - 1)M 2 

Z CO 


HI 


h/(h-1) 


Reduced pressure coefficient, = 


£(h+ 1)M 2 ] 1/3 


77 * 


- 1 


_ 1 2 

Pressure ratio, - « h M c +1 

2 00 p 


Machnumber, -M 


co u 



r / \ 2 1 


i + * : ' m 2 

i-^\ 


2 ® 

\ U / 



L \ »/ J 



- 1/2 


Critical Machnumber, - M* 


The critical value of the pressure coefficient (C*) is evaluated 
from the exact relation. 


Note: 

CP(STAR) 
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FIGURE 3 - Geometric division of an arbitrary planform into a number of 
trapezoidal panels. 
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FIGURE 6 - Block diagram of the computer program. 
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** Subcritical flow past a lifting airfoil. 
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FIGURE 12 - Discontinuous flow past a non-lifting airfoil. 
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FIGURE 15 - Comparisons of analytical and numerical results for a symmetrical wing. 
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FIGURE 16 - Incompressible flow solutions for high aspect ratio wings 


with a blunt leading edge. 
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FIGURE 17 - Subcritical flow past a non-lifting rectangular high aspect ratio wing. 
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FIGURE 22 - Discontinuous flow past a non-lifting rectangular wing 
near the mid-span station. 
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FIGURE 24 - Supercritical flow past a lifting rectangular wing near the mid-span station 
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FIGURE 26 - Flow chart for the calling sequence of the program subroutines. 
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